NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


THESIS 

A  NUMERICAL  STUDY  OF  HEAT  TRANSFER  BEHAVIOR 

IN  WELDING 

by 

Yasar  Vehbi  Isiklar 
June,  1998 

Thesis  Advisor:  Ashok  Gopinath 

Approved  for  public  release;  distribution  is  unlimited. 


[SffIC  QUALITY  K^CPECTBD  1 


19980803  032 


REPORT  DOCUMENTATION  PAGE 

Form  Approved  0MB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215 

Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperworic  Reduction  Project  (0704-0188)  Washington  DC 

20503. 

1.  AGENCY  USE  ONLY  (leave  Wa«ifc;  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERH) 

June  1998  Master’s  Thesis 

4.  TITLE  AND  SUBTITLE 

A  NUMERICAL  STUDY  OF  HEAT  TRANSFER  BEHAVIOR  IN  WELDING 

5.  FUNDING  NUMBERS 

6.  AUTHOR  (S)  Isiklar,  Yasar  V. 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

8.  PERFORMING 

ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

NAVSEA  Carderock,  Maryland 

10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

1 1 .  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the 
official  policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 

12a.  DISTREBUnON/AVAILABlLITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (maximum  200  words) 

A  numerical  model  has  been  developed  for  three-dimensional  transient  conduction  based  temperature  calculations  in 
underwater  wet  welding  on  a  thick  rectangular  plate.  The  numerical  scheme  is  based  on  a  fully  implicit  finite  volume 
method.  A  variable  mesh  size  centered  around  the  moving  heat  source,  and  temperature  dependent  thermal  properties  have 
been  used  in  the  calculations.  Convective,  radiative  and  boiling  surface  thermal  conditions  have  also  been  included.  The 
weld  pool  region  itself  has  been  modeled  as  a  solid  region  of  thermal  conductivity  higher  than  the  surrounding  unmelted 
region.  The  validity  of  the  results  was  checked  by  comparison  with  Rosenthal's  three-dimensional  solution  for  a  moving 
point  heat  source,  and  other  results  in  the  literature. 

14.  SUBJECT  TERMS  Underwater  Wet  Welding,  Heat  Transfer,  Finite  Volume  Method 

15.  NUMBER  OF 

PAGES  125 

16.  PRICE  CODE 

17.  SECURITY  CLASSmCA- 
TTON  OF  REPORT 
Unclassified 

18.  SECURITY  CLASSIFI¬ 
CATION  OF  THIS  PAGE 
Unclassified 

19.  SECURITY  CLASSIFICA¬ 
TION  OF  ABSTRACT 
Unclassified 

20.  LIMITATION  OF 
ABSTRACT 

UL 

NSN  7540-01-280-5500  Standard  Form  298  (Rev.  2-89) 


Prescribed  by  ANSI  Std.  239-18  298-102 


i 


11 


Approved  for  public  release;  distribution  is  unlimited. 


A  NUMERICAL  STUDY  OF  HEAT  TRANSFER  BEHAVIOR  IN  WELDING 


Yasar  V.  Isiklar 

Lieutenant  Junior  Grade,  Turkish  Navy 
B.S.,  Turkish  Naval  Academy,  1992 

Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  MECHANICAL  ENGINEERING 

from  the 


NAVAL  POSTGRADUATE  SCHOOL 
June  1998 


Author: 


Approved  by: 


iv 


ABSTRACT 


A  numerical  model  has  been  developed  for  three-dimensional  transient 
conduction  based  temperature  calculations  in  underwater  wet  welding  on  a  thick 
rectangular  plate.  The  numerical  scheme  is  based  on  a  fully  implicit  finite  voliune 
method.  A  variable  mesh  size  centered  arormd  the  moving  heat  source,  and  temperature 
dependent  thermal  properties  have  been  used  in  the  calculations.  Convective,  radiative 
and  boiling  surface  thermal  conditions  have  also  been  included.  The  weld  pool  region 
itself  has  been  modeled  as  a  solid  region  of  thermal  conductivity  higher  than  the 
surrounding  unmelted  solid  region.  The  validity  of  the  results  was  checked  by  comparison 
with  Rosenthal's  three-dimensional  solution  for  a  moving  point  heat  source,  and  other 
results  in  the  literature. 
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I.  INTRODUCTION 


To  improve  the  quality  of  the  underwater  welding  and  to  accomplish  a  reliable, 
permanent  underwater  wet  welding  capability  has  a  great  importance  in  today’s  industrial 
and  military  facilities.  With  the  development  of  underwater  wet  welding  techniques,  the 
time  and  the  money  required  for  permanent  and  temporary  repairs  of  ships  and  other 
underwater  structures  can  be  minimized.  Today,  the  use  of  hyperbaric  welding  process 
obtains  a  limited  quality  of  welding  especially  for  the  construction  and  repairs  of 
underwater  pipelines.  In  this  process  a  large  pressure  chamber  is  used  to  keep  the  water 
away  from  the  workpiece.  But,  operating  this  kind  of  chamber  is  very  expensive  and  due 
to  the  limited  geometric  size,  only  a  few  joint  configurations  can  be  enclosed  in  a 
chamber  [Refl].  The  other  welding  techniques  such  as  double  shielding  and  flux 
shielding  also  use  the  water  removing  theory  from  the  arc  area  during  welding.  But,  the 
working  area  must  be  completely  prevented  from  water  for  satisfactory  welding  results 
[Ref  2]. 

Currently,  underwater  wet  welding  is  used  for  the  temporary  repair  needs. 
Because  of  their  poor  quality  compared  to  surface  (air)  welds  (they  obtain  80%  of  the 
tensile  strength  and  50%  ductility  of  the  surface  welds  [Ref  3].),  they  must  be  replaced  as 
soon  as  possible.  Therefore,  the  development  of  a  more  efficient  wet  welding  technique  is 
the  only  solution  to  this  problem. 

The  surrounding  water  environment  during  wet  welding  causes  rapid  cooling  and 
steep  temperature  gradients  in  the  weld  area  behind  the  arc  [Ref. 4].  Because  of  the 
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extremely  complex  nature  of  the  heat  transfer  phenomena  between  heated  surface  and  the 
water  environment,  a  numerical  model  simulation  is  necessary. 

In  the  present  study,  a  numerical  model  has  been  developed  for  transient,  three- 
dimensional  conduction  heat  transfer  in  imderwater  welding  process  on  a  thick 
rectangular  plate.  The  numerical  scheme  was  based  on  fully  implicit  finite  volume  model, 
including  convection,  radiation  and  boiling  surface  thermal  boundary  conditions.  The 
different  regimes  of  boiling  were  accounted  for  on  the  surface.  A  variable  mesh  size 
centered  around  an  arc  source  moving  at  constant  speed  was  used  to  determine 
temperature  variations  inside  and  aroimd  the  weld  pool.  The  weld  pool  region  itself  has 
been  modeled  as  a  solid  region  of  thermal  conductivity  higher  than  the  surrounding 
unmelted  solid  region.  The  input  data  from  the  previous  studies  based  on  different 
methods  were  used  to  check  the  accuracy  and  the  validity  of  the  numerical  method. 


2 


II.  BACKGROUND 


A.  PREVIOUS  STUDIES 

Rosenthal  did  the  most  important  early  work  on  the  theory  of  the  effect  of  moving 
sources  of  heat  in  the  late  of  1930s.  He  studied  the  fundamentals  of  this  theory  and 
derived  equations  for  two-dimensional  and  three-dimensional  heat  conduction  in  a  solid 
when  a  moving  source  is  in  use  [Ref.  5,6].  The  analytical  solution  derived  by  Rosenthal 
was  based  on  the  principle  of  a  quasi-stationary  thermal  state.  The  quasi-stationaiy 
thermal  state  represents  a  steady  thermal  response  of  the  weldment  with  respect  to  the 
moving  coordinates.  In  other  words,  the  origin  moves  with  the  heat  source,  and  for  an 
observer  at  origin,  the  temperature  distribution  and  the  pool  geometry  do  not  change  with 
time  [Ref.  6,7].  Rosenthal's  solution  for  three-dimensional  heat  flow  during  welding  is 
as  follows  [Ref.  6]; 


l7iiT-T,)k^R 

Q 


=  exp 


^  U{R-x)^ 

2a, 


(2.1) 


where 

(  2  2  2V^2 

X  +}>  +z  ) 

r  =  temperature, 

Tq  =  workpiece  temperature  before  welding, 

=  thermal  conductivity  of  solid, 

Q  =  heat  input  to  the  workpiece. 
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U  =  welding  speed, 

X  =  the  distance  that  the  heat  source  traveled, 

a,  =  thermal  diffusivity  of  solid  (i.e.,  /  pC^ ,  where  p  and  are  density  and 

specific  heat  of  solid,  respectively). 


The  assumptions  used  by  Rosenthal  to  derive  the  equation  above  are  as  follows 
[Ref  7]: 

1 .  Point  heat  source. 

2.  No  melting  and  negligible  heat  of  fusion. 

3.  Constant  thermal  properties. 

4.  No  heat  loss  from  the  workpiece  surface. 

5.  Infinitely  wide  workpiece. 

Although  Rosenthal's  assumptions  help  to  simplify  the  mathematical  analysis 
involved,  there  are  some  significant  deviations  between  theoretical  and  experimental 
results.  For  instance,  as  a  result  of  the  point  heat  source  assumption,  the  temperature  at 
the  weld  centerline  goes  to  infinity  even  though  the  power  of  the  heat  source  is  finite. 
Also,  the  values  of  thermal  properties  change  with  the  temperature  and  neglecting  the 
heat  fusion  gives  considerable  errors.  [Ref  7] 

Tanaka  did  the  first  studies  on  the  application  of  the  mathematical  analysis  of 
non-stationary  heat  flow  to  the  practical  problems.  Naka  and  Masubuchi  studied  the 
mathematical  analysis  of  non-stationary  heat  flow.  They  used  dimensional  expressions  to 
make  the  numerical  calculations  simple.  Nippes  and  Savage  studied  the  cooling  rates  of 
heat  affected  zones  by  using  a  graphical  approach  method.  Suzuki  found  an  analytical- 
empirical  method  of  studying  the  effects  of  welding  parameters  and  determined  the 
cooling  rate  from  welding  conditions  with  the  help  of  a  monograph.  [Ref.6] 
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Adams  derived  new  equations  from  Rosenthal’s  equation  by  using  the  fusion  line 
as  a  boundary  condition  and  calculated  the  peak  temperature  at  a  distance  from  the  fusion 
boundary  at  the  weldment  surface  [Ref.  7].  His  equation  for  three-dimensional  heat  flow 
is  as  follows: 


1  _  SAAnk^a^ 

Tp-T,~  QU 

where  all  the  terms  are  defined  in  equation  (2.1)  except 

Tp  =  peak  temperature  at  a  distance  Y  from  the  fusion  boundary, 

=  melting  temperature, 

Christensen  also  used  analytical-empirical  approach  and  derived  dimensionless 
equations  based  on  Rosenthal's  three-dimensional  equation.  He  explained  the  relation¬ 
ships  between  the  welding  conditions  and  the  weld  bead  geometry.  [Ref.  7] 

Recently,  numerical  analysis  methods  and  computer  programs  have  been 
commonly  used  to  develop  the  previous  assumptions.  In  1965,  The  Battelle  Institute 
Geneva  Laboratory  in  Switzerland  conducted  a  computer-aided  study  about  analyzing  of 
heat  flow  in  weldments  [Ref  6].  The  University  of  Wisconsin  and  McDonnell  Douglas 
Aircraft  Company  also  conducted  similar  studies  in  the  same  year  [Ref.  6].  In  1970, 
M.I.T.  researchers  studied  heat  flow  during  underwater  welding.  They  also  developed 
finite  element  programs  on  heat  flow  during  welding  [Ref  6].  Oreper  and  Szekeley 
examined  stationary,  axisymmetric  TIG  (tungsten-inert-gas)  welding  process  with  a 


2  + 


\^^sj 


+  - 


Tm-To 


(2.2) 
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moving  boundary  by  using  finite  difference  method.  Their  formulation  contained  the 
affects  of  transient  conduction,  electromagnetic,  buoyancy  and  surface  tension  forces 
[Ref  8].  Kou  and  Wang  performed  computational  studies  of  the  GTA  welding  process. 
They  presented  a  computer  simulation  of  three-dimensional  convection  for  an  arc  source 
moving  at  constant  speed.  They  considered  electromagnetic,  buoyancy  and  surface 
tension  forces  on  the  pool  surface.  They  found  very  good  agreement  between  the 
calculated  and  observed  fusion  boundaries  [Ref  9].  They  also  studied  the  computer 
simulation  of  three-dimensional  convection  in  laser  melting  by  considering  the  buoyancy 
force  and  the  surface  tension  gradient  at  the  weld  pool  surface  [Ref.  10].  Correa  and 
Sundell  studied  axisynunetric  stationary  arc  source  by  using  different  grid  sizes  for 
computation  of  flow  and  electromagnetic  fields  [Ref  11].  Saedi  and  Unkel  developed  a 
thermal-fluid  model  of  the  weld  pool.  Their  model  was  based  on  the  stationary  arc.  To 
describe  the  weld  pool  geometry,  they  matched  the  convective  and  the  conductive  heat 
fluxes  at  the  weld  boundary  by  using  an  iterative  calculation  method  [Ref.  12].  Zacharia 
et  al.  made  three-dimensional  calculations  on  the  effects  of  the  heat  source  in  the 
stationary  GTAW  process.  He  indicated  that  a  depressed  area  formed  at  the  weld  pool 
center  because  of  an  outward  fluid  flow  caused  by  surface  tension  force  [Ref  13].  Kim 
and  Na  developed  a  model  on  heat  and  mass  flow  for  stationary,  GTAW  process  with 
electromagnetic,  buoyancy  and  surface  tension  forces.  They  used  numerical  mapping 
method  for  calculations  [Ref  14].  Ramanan  and  Korpela  compared  the  effects  of  thermo¬ 
capillary  and  Lorentz  forces  on  the  flow  pattern  in  a  stationary  weld  pool  with  buoyancy 
forces.  They  used  multi-grid  methods  and  a  local  grid  refinement  technique  with 
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axisyimnetric  stationary  arc  source  [Ref.  14].  Ule,  Joshi  and  Sedy  determined  three- 
dimensional  transient  temperature  variations  in  the  GTAW  process  hy  using  the  explicit 
finite  difference  method.  They  used  different  mesh  sizes  and  temperature  dependent 
thermal  properties.  They  also  considered  convective  and  radiative  surface  thermal 
conditions  during  calculations  [Ref.  16].  Kanouff  and  Greif  studied  the  unsteady 
development  of  an  axisymmetric  arc  weld  pool  in  GTAW  process.  They  used  moving 
grids  to  follow  the  phase  change  boundary  and  considered  the  effects  of  Marangoni, 
Lorentz  and  buoyancy  forces  in  the  calculations  [Ref.  17].  Joshi,  Dutta,  Schupp  and 
Espinosa  developed  a  three-dimensional  numerical  model  to  describe  the  flow  circulation 
phenomena  in  aluminum  weld  pools  under  non-axisymmetric  Lorentz  force  field  [Ref. 
18,19]. 


B.  FINITE  VOLUME  METHOD 

The  finite  volume  method  is  one  of  the  simple  and  well-established 
Computational  Fluid  Dynamics  (CFD)  techniques  that  were  originally  developed  as  a 
special  finite  difference  method  [Ref.  19].  The  stages  of  the  numerical  algorithm  in  this 
method  are  as  follows  [Ref.  19]: 

1.  Formal  integration  of  the  governing  equations  of  fluid  flow  over  all  the  finite 
control  volumes  of  the  solution  domain. 

2.  Discretisation  involves  the  substitution  of  a  variety  of  finite  difference  type 
approximations  for  the  terms  in  the  integrated  equation  representing  flow 
process  such  as  convection,  diffusion  and  sources.  This  converts  the  integral 
equations  into  a  system  of  algebraic  equations. 

3.  Solution  of  the  algebraic  equations  by  an  iterative  method. 
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In  the  control  volume  integration,  the  calculation  domain  is  divided  into  discrete 
control  volumes.  There  is  only  one  control  volume  surrounding  each  grid  point  and  the 
boundaries  of  the  control  volumes  are  positioned  side  to  side  in  the  middle  of  the  distance 
between  the  grid  points.  Integrating  the  governing  differential  equation  over  each  control 
volume  derives  the  resulting  discretised  equations.  The  discretised  equations  express  the 
conservation  of  quantities  such  as  mass,  momentum  and  energy  for  each  control  volume. 
This  characteristic  exists  for  any  number  of  grid  points.  The  numerical  solution  insures 
the  validity  of  the  conservation  principle  over  the  whole  calculation  domain  for  the 
related  quantities.  Versteeg  and  Malalasekera  [Ref  20]  wrote  "  This  clear  relationship 
between  the  numerical  algorithm  and  the  underlying  physical  conservation  principle 
forms  makes  finite  volume  method  much  simpler  to  imderstand  by  engineers  than  finite 
element  and  spectral  methods".  To  understand  the  finite  volume  method  better  an 
illustrative  example  can  be  given  for  one-dimensional  steady  state  heat  conduction 
situation.  [Ref  20,21,22] 

1.  One  dimensional  steady  state  heat  conduction 

The  governing  equation  for  one-dimensional  steady  state  heat  conduction  is 


dx\  dx ) 


+  5  =  0 


where 

k  =  thermal  conductivity, 

7’=  temperature, 

S  =  the  rate  of  heat  generation  per  unit  volume. 


(2.3) 
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The  domain  is  divided  into  small  and  nonoverlapping  control  volumes.  A  part  of 
one-dimensional  grids  generated  is  shown  in  Figure  2.1.  Here,  The  grid  point  under 
construction  is  denoted  as  P  and  the  neighboring  nodes  to  the  east  and  west  are  denoted  as 
E  and  W  respectively.  The  lower  case  letters  e  and  w  denotes  the  east  and  west  faces  of 
the  control  volume.  The  distance  between  the  nodes  W  and  P  is  denoted  by  and 

between  P  and  E  is  denoted  by  Scp^ .  The  distances  between  w  and  P  and  between  P  and  e 
are  given  by  Sc„p  and  Scp^  respectively.  The  control  volume  width  is  shown  as 
Ax  =  (S)c^, 


Figure  2.1  One-dimensional  grid  [Ref.  20] 


Integrating  equation  (2.3)  over  the  control  volume  gives. 


fr 

k — 


(2.4) 
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By  assuming  a  piecewise-Iinear  profile  assumption  (Figure  2.2),  equation  (2.4) 


can  be  evaluated  as 


Te-Tp 


6x 


-k.. 


f  T  —T  ^ 

P 


PE  J 


+  5Ax  =  0 


WP  J 


(2.5) 


where  S  is  the  average  value  of  S  over  the  control  volume. 


Figure  2.2  Piecewise  -linear  Profile  [Ref.  22] 


By  arranging  the  terms,  the  final  discretised  equation  can  be  written  as 

OpTp  =  cipT^  +  cifyT^  +  b 


where 


a 


E  ~ 


dXpE 


(2.6) 


(2.7a) 
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(2.7b) 


aj,=a^+  Ofy 


b  =  SAx 


(2.7c) 


(2.7d) 
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in.  MODEL  DEVELOPMENT 


A.  DEFINING  THE  WORKPIECE 

The  workpiece  has  been  defined  as  a  thick  rectangular  plate  (Figure  3.1). 


Figure  3.1  The  workpiece 


where 

1 .  Right  Lateral  Face  (East) 

2.  Left  Lateral  Face  (West) 

3.  Back  Face  (North) 

4.  Front  Face  (South) 

5.  Top  Face  (Top) 

6.  Bottom  Face  (Bottom) 

Convective  heat  transfer  coefficients  for  each  face  have  been  defined  as  follows: 

Right-Lateral  Face: 

Left-Lateral  Face  : 

Back  Face  :  h„ 

Front  Face  : 

Top  Face  :h, 

Bottom  Face  :  hf, 

The  same  method  has  been  also  used  to  define  thermal  conductivity  and  heat  flux  terms 
for  each  face. 
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B.  BOUNDARY  CONDITIONS 


1.  Top  Face 

By  adding  radiative  and  convective  heat  loss  from  the  top  face,  boundary 
condition  equation  can  be  defined  as 

=  -r.)+?"+<74rl  -c)  O.i) 

In  equation  (3.1),  for  the  radiation  term,  can  be  neglected.  »  C).The 

distance  between  the  node  point  P  and  the  face  can  be  taken  as  Ax/2,  where  the  distance 
between  the  node  point  P  and  the  neighboring  nodes  is  Ax .  Again,  for  the  radiation  term, 
because  of  the  small  distance  between  the  node  point  P  and  the  face,  by  assuming 
s  Tp ,  we  have 


-k 


inwall  ) 
lSxl2 


=  KT^an  -TJ  +  q"+CT£T^ 


(3.2) 


Equation  (3.2)  can  be  opened  as 


Ax 


^  wall 


=  hT^.^„-hT^+q”+CTsT; 


(3.3) 


Equation  (3.3)  can  be  arranged  as 
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(3.4) 


^wall 


(2k 
—  +  h 


--T,+hT„-q"+asT; 

Ax 


Equation  can  be  written  as 


^  wall 


Ifc 

—T,+hT^-q"+asT^ 

Ax 


2k 
—  + 
Ax 


h 


(3.5) 


The  unknown  temperature  for  the  top  face  points  is  found  as 


^wa/l 


2kTp  +  hAxT„  -q" Ax- asAxT* 
2k  +  hAx 


(3.6) 


where 


<l'\cp  =  ^"source+^",+<l'\ciImg 


(3.7) 


q" source  souTce  heat  flux, 

q'\  =  a  constant  arbitrary  heat  flux  which  may  be  applied  to  the  top  face, 

=  boiling  heat  flux. 


2.  Other  Faces 

By  using  the  same  method  from  the  top  face  (without  radiation). 


15 


(3.8) 


2kTp  +h^xT^  -q"Ax 
2k  +  /zAx 


3.  Simulating  the  Arc 

To  simulate  the  heat  input  from  the  arc  to  the  workpiece,  it  is  assumed  that  the 
heat  input  distribution  of  the  arc  have  a  Gaussian  distribution  on  the  top  face  of  the 
workpiece.  The  general  equation  is 

“  -4''^ 

Q  =  '^^dr  (3.9) 

0 

where 


Q  =  the  total  heat  input  into  the  workpiece, 

^0  =  the  volumetric  energy  generation  rate, 

^0  =  the  radius  of  the  heat  input  distribution, 
d  exponential  factor, 

By  solving  equation  (3.9),  the  volumetric  energy  generation  rate  can  be  found  as 


— 


Qd 


TO-. 


(3.10) 


4.  Boiling  Heat  Transfer 

Modes  or  regimes  of  boiling  and  the  related  equations  can  be  classified  as  follows 
(where  AT,  =7;- r,^,): 
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a. 


Free  Convection  Regime  (Ar^  <  5  °C) 


In  this  regime,  natural  convection  effects  determine  the  heat  transfer 
between  the  heating  surface  and  surrounding  liquid.  Recommended  correlations  for  upper 
surface  of  heated  plate  are  [Ref.  25] 


Nul  =  0.54P<' 

(lO^  <10’) 

(3.11) 

IhiL  =0.15Ral'^ 

(lO’  <Ra^^  <10") 

(3.12) 

where  the  Rayleigh  number , 


va 


(3.13) 


here 


g  =  gravitational  acceleration,  m/s 
P  =  volumetric  thermal  expansion  coefficient,  K'^ 

V  =  kinematic  viscosity,  m  /s 
a  =  thermal  diffiisivity,  m^/s 

L  =  characteristic  length,  L  =  Plate  surface  area  {A^)l  Perimeter  (P) 


7  Nuik 

h  = - 

L 


(3.14) 


and  the  value  of  heat  flux  is. 
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q"=KT,-Tj 


(3.15) 


b.  Nucleate  Boiling  Regime  (5  °C  <  AT^  <  30  "C) 

The  most  useful  nucleate  pool  boiling  correlating  equation  was  developed 
by  Rohsenow  [Ref.  23,  24,  25,], 


qs'=Mihfg 


SiPl  -Pv)' 

1/2 

‘.AT,  ] 

(T 

(3.16) 


where  the  subscripts  s,  1,  and  v  express  surface,  saturated  liquid  state  and  vapor  state.  The 
definition  of  each  term  in  equation  (3.16)  are  as  follows: 


Hi  =  viscosity  of  the  liquid,  kg/ms 
=  latent  heat  of  vaporization,  J/kg 
g  =  gravitational  acceleration,  m/s^ 

=  density  of  the  saturated  liquid,  kg/m^ 

=  density  of  the  saturated  vapor,  kg/m^ 

^  =  surface  tension  of  the  liquid-to-vapor  interface,  N/m 

Cp  ,  =  specific  heat  of  saturated  liquid,  J/kgK 

AT,  =  T, 

Pr,  =  Prandtl  number  of  the  saturated  liquid 
n  =1.0  for  water,  1 .7  for  other  fluids 

C^  j  =  empirical  constant  that  depends  on  the  nature  of  the  heating 

surface  fluid  combination  and  whose  numerical  value  varies 
from  system  to  system 


But,  in  underwater  welding,  the  surrounding  water  temperature  is  below  the 
saturation  temperature  (between  0  °C  and  30  °C ).  This  is  called  as  the  heat  transfer  to  a 
subcooled  liquid.  For  the  subcooled  boiling,  the  heat  flux  can  be  estimated  as  [Ref  23] 
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(3.17) 


r 

.  1  4. 

7^/  ^sat  '^liquid  ) 

24 

r  2  1 

Pv 

1/4" 

^TraiT 

_Og{pl-py)_ 

where 

k 

‘  =  thermal  conductivity  of  the  liquid,  W/m.K 

=  thermal  diffusivity  of  the  liquid  (A://x;^),m^/s  and. 


T  = 


a 

1/2 

2  "I 

Pv 

.g(Pl  -Pv). 

_Og(pl-py)_ 

(3.18) 


c.  Transition  Boiling  Regime  (30°C  <  AT^  <  120  "C) 

For  the  transition-boiling  regime,  no  sufficient  theory  has  been  derived. 
This  regime  is  between  the  maximum  and  minimum  heat  fluxes  where  [Ref.  23,25], 


?"max  =  0-149/j^gA 


ogU-Pv) 


ni/4 


(3.19) 


^"nun  =  0.09p^hj.g 


ga{p,-p^) 

.  iPi+P.f 


nl/4 


(3.20) 


By  assuming  a  linear  heat  flux  distribution  in  the  transition  boiling  regime 
(Figure  3.2),  it  can  be  written. 


l0g(g"max  )-10g(^"min  )  l0g(g")  "  l0g(g"min  ) 

Iog30-logl20  log ATg -log  120 


(3.21) 
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by  arranging  equation  (3.21), 


log 


log(9")  =  - 


^  max 

f  'i.C\\ 


AT 


log 


30 

Vl20y 


120 


log  T-r  +\og{q"^J 


(3.22) 


and  the  heat  flux  at  a  point  in  transition  boiling  regime  can  be  found  as 


^”=10'°®^^"^ 


(3.23) 


Figure  3.2  Linear  Distribution  of  Heat  Flux  in  Transition  Regime 

d.  Film  Boiling  Regime  (AZg  >  120°  C) 

In  film  boiling  regime  for  flat  horizontal  surfaces,  Westwater  and  Breen 
recommended  the  following  correlation  for  conduction  heat  transfer  coefficient  [Ref  23], 
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he  =  0.59< 


giPi  -p^)pX^\hfs  +0-68c^,ArJ] 


^PXTe 


1/4 


where 

=  viscosity  of  the  vapor,  kg/ms 
=  thermal  conductivity  of  vapor,  W/mK 
Cp,  =  specific  heat  of  saturated  vapor,  J/kgKand, 


X  =  2n\ 


G 


,1/2 


g{Pl-P.)\ 


(3.24) 


(3.25) 


Bromley  suggested  combining  conduction  and  radiation  heat  transfer  coefficients 
[Ref  23], 


htotal  — hc'i'O.lShr 


where 


rp4  _  rp4  \ 
s  ^  sat 


here 

* 

=  surface  emissivity 
=  absolute  surface  temperature 

and  the  resulting  heat  flux  in  this  regime  is  found  as 


^"=  htotal 


(3;26) 


(3.27) 


(3.28) 
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C.  COEFFICIENTS  USED  IN  THE  EQUATIONS 


The  coefficients  used  deriving  the  discretised  equations  are  as  follows: 


dx 


PE 


Qfy  — 


dx. 


WP 


= 


k„A„ 

n  n 


PN 


M 

^SP 


a  =  ^  ^ 


Or  = 


KA 


^  dz 


FT 


^BP 


0  AV 


Op  =  &{ap  +a^  +  Qp!  +  +  flj.  +  ap)+  a 


coeffe= 


2k ^  4-  h^Sxp^ 


coeffiv= 


2k^.+h^.dxjyp 


coeffii  = 


^n^PN 


^k„+K^p^, 


coeffi  = 


2k ^  +h,dzpp 


fliacoeffe  = 


dx 


PE 


2k^  +  hg&Cpp 


coeffs  = 


K^sp 


2k  ^  +h^dysp 


coefflT= 


2kf^  +  hj^dz  pp 


flvxcoejfw  - 


dx 


WP 


2k  ^  +  h^dXfyp 
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fluxcoejfn 


^PN 


2k „  + 


fluxcoejft  = 


dz 


PT 


2k,  +h,dzpj. 


radcoeff  = 


2k,  +h,Szpj. 


fluxcoeffs  = 


^SP 


2k ^  +h^dysp 


flvxcoeffb  = 


dz 


BP 


2k  I,  +  h,,dzpp 


D.  DERIVATION  OF  THE  EQUATIONS 

The  problem  is  governed  by  the  equation 


cbc ) 


Qy  j 


d  (,  dT\ 
+  —  k — 


1.  Interior  Nodes 


(3.29) 
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Integration  of  equation  (3.29)  over  the  control  volume  and  over  a  time  interval  from  t  to 
t  +  Ar  gives 


(3.30) 


This  may  be  written  as 


Adxdt-k- 


t+M 


1 1-  *- 


Adydt 


/+A/ 


/  cv 


dz  ) 


Adzdt 


(3.31) 


where  A  is  the  face  area  of  the  control  volume  and  dx,  dy,  dz  are  the  dimensions  of  the 
control  volume.  Equation  (3.31)  may  be  written  as 


i 


'7  ' 

pc — dt 
J  ^  dt 


CH.  / 


dv=  I 


t  L 


kA 


dx 


4^ 

I  a. 


r+A/ 

dt+  I 


kA 


kA 


dt 


/+A/ 

r 

f, 

( ,  .57) 

1 

kA  — 

kA — 

j 

t 

A  ) 

t 

1  5z  j 

dt 


(3.32) 
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By  assuming  the  grid  point  value  of  the  temperature  at  a  node  to  prevail  over  the  whole 
control  volume,  the  left  hand  side  of  equation  (3 .32)  can  be  written  as 


(3.33) 


where  Tp  refers  to  temperatures  at  time  t  and  Tp  refers  to  temperatures  at  time  t  +  M. 

By  applying  central  differencing  to  the  diffusion  terms  on  the  right  hand  side  equation 
(3.32)  can  be  written  as 


t+6t 

ixi^,-nw=  \ 


f  yr 


Sc, 


( 


•PE  J 


n-Tu 


KA  - 


w 


r-fAr 

cit+  J 


f  7 

U  " 


FN  J 


44 


T  -T'\ 
fp _ fs 

^SP  J 


^P 

&PT  J 


44 


Tp-T, 


'BP 


(3.34) 


The  values  of  Tp,Tp,T„,,Tj^,T^,Tp  and  Tg  vary  with  time.  The  time  integral  can  be 

calculated  by  using  temperatures  at  time  t  or  at  time  r  +  Ar  or,  a  combination  of 
temperatures  at  time  t  and  r  +  Ar.  This  approach  may  be  generalized  by  defining  a 
weighting  parameter  ©  between  0  and  1 . 


rf(=[er, +(i-0)r,°]i» 


t 


(3.35) 
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In  equation  (3.35)  for  ©  =0  the  temperature  at  time  level  t  is  used  (Explicit  scheme);  for 
©  =0.5  the  temperatures  at  t  and  /  +  Arare  equally  weighted  (Crank-Nicholson  scheme); 
for  ©  =1  the  temperature  at  time  level  r  +  A/  is  used  (Fully  implicit  scheme).  By  using 
equation  (3.35)  and  dividing  A/  throughout,  equation  (3.34)  may  be  arranged  as 


Ar 


M,  ^ 

^PE  ^WP  ^PN  ^SP 


+  k.A. 


An-Ts) 


& 


-k.A 


b^b 


PT 


dz 


BP 


+  (1-©) 


Sx 


PE 


^PN  ^SP 


dz 


PT 


dz 


BP 


(3.36) 


Equation  (3.36)  may  be  re-arranged  as 


fc—+e  Tp=^[qt^ 

L  A?  ^^P  ^/>yV  ^SP  ^PT  ^BP  J  ^PE 


+(i-0K]+^K  +(i-©W1+^K +{i-0K]+|4-[©r, +(i-0)7j] 

^^P  ^PN  ^SP 

+(\-®xV^\^, +(i  -0K] 


& 


'PT 


BP 


AF  ( 

pc - (1-0) 

^  Ar  ^ 


^PE 


kA  kA^  k.A,  krA. 


^PN  ^SP  ^PT  ^ 


BP  J 


i  p 


(3.37) 
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By  putting  the  defined  coefficients  from  the  previous  section,  equation  (3.37)  may  be 
written  as 

[a®  +©(«£  +a^  +as  +a^  +  )]?>  =  +(l-©)r£] 

+ [©r^  +  (1  -  ©)r;  ]+  a,  [©7;  +  (1  -  ©W  ]+ a,  [©r,  +  (1  -  ©fc®  ] 

+ a,  [©r, + (1  -  ©)r,®  ]+ a  ^©T’g  +  (1  -  ©)r;  ] 

+  [a®  -  (1  -  ©)(a£  +a^  +af^  +as+a.j.+ag  )]rp®  (3.38) 


Finally,  by  grouping  the  known  and  the  unknown  terms  at  each  side,  the  discretised 
equation  is  foimd  as 


apTp  -®apTp  =  (1-©) 

[a,T^  +a^T^  +aX  +^s^s  +apT^]+[al  -(l-0)(t^£  +  % 

+  +  ag  + Cp  +  Cgi^p  (3.39) 


2.  Left-Front-Top  Comer 
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By  using  the  same  method  from  equations  (3.30-3.33),  equation  (3.29)  can  be  written  as 


pdj-,-ltW=  ] 


T  -T 


5c 


'PE  J 


r+A/ 

7  Tr-TA 

f  T  -T  ^ 

+  f 

k,A  ^  ^ 

KA 

J 

t 

.1  &PT  ) 

0  D 

^  ^BP  J_ 

•~r  _ fr  \  f+iV 

^  ! 

V  J  t 


\it 


KA 


T,-T, 

^PN  . 


^SP  / 

(3.40) 


[it 


In  equation  (3.40),  and  take  the  value  of  equation  (3.8)  and  Tj  takes  the  value  of 
equation  (3.6).  By  using  these  values  and  equation  (3.35),  we  have 


fX 


A/ 


Sx 
2k, A 


(r, 


^PN  ^SP 


PE 

r 


&c, 


WP 


Tp- 


2k  Jp  +  KSx^T„  -  q\ 


2k  +  h^&c^ 


Tp- 


2k  Jp  +A^spT^  -q'\  dysp 


SP 


2k, A 


dz 


PT 


'Ik^Tp +h^& ppT^  q  iQpSzpj.  crs&ppTf 

2k^  A-h^dZpj 


-T 


KA 


dz 


Jp-Tp) 


BP 


+  (1-0)1 

k^A 


^kJTp  h^Sx^rpT^  Q  w^\ 


WP 


Sk 


PE 


Sx 


WP 


2k^,  +  hSx, 


2k  A 


r 


Sy 


PN 

2kA 


SP 


w^WP 


I  »  »  I  j-0  _  ^^s'^P  K^Sp'^o:>  ^  s  ^SP 


2k,  +h,dy^ 


SP 


dz 


PT 


2k Jl  +h,&pj.T^  -q'\^p  -(^s&ptT 
2k,  +h,dziyj- 


T'O 
p 


J 

-Tl) 


dz 


BP 


(3.41) 


Equation  (3.41)  may  be  arranged  as 
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2k A  ( K^spT. -q"s  ^sA  [  ^^prT. 

^SP  V  "^^s^SP  J  &PT  i 


-(j'top&PT-^^S&pTTp 

2k,  +h,dzpp 


(3.42) 


Equation  (3.42)  may  be  re-arranged  as 


I  I  2^.4.r  l>A^  V.4 , 2A.4[  '' 

At  _^PE  ^WP  \^^w'^\^wpj  ^PN  ^SP  \^K’^K^SP  J 


KS>Cw^^  q'\s>Cfpp  \^2kA(  <i",^sp  ^ ,  2M 

2k„  +  h,,&jpp  2k^  +  h„Sc^  J  ^sp[.'^s  +  K^SP  +  K^SP  J  ^PT 
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h,&pj.T^ _ q  top  &PJ-  (xe&pjTp  ^  ^  ^  0n  A 

2,k^  +hi&pp  2kj  +h^&pp  2k ^  +h,&pp  j  Sxy^p  '^^w^ivp 

q\dK^  Y2k^A(  h^^,pT^  q\Sy,p  ]  ,  2^,4  T  h,dzppT^ 

'2.k^  J  ^ sp  Kp'^s  "^^s^SP  "^^s^SP  J  ^PT  +h^&pp 


top  ^PT  creSzppTp  ^ 

2k,  +  h,Szpp  2k,  +  h,dzpj  J 


By  putting  the  defined  coefficients  from  the  previous  section,  equation  (3.43)  may  be 
arranged  as 


[al  +  ©[a^  +  2a^  (coe^)  +  a^+  20^  icoeffs)  +  2aj  {coeffi)  +  jjfp 

=  [eTp  +  (1  -  ©)r;  ] + [qt,  +  (1  -  e)T^  ] + [©r,  +  (i  -  @)t^  ] 

+  [a°p  -  (1  -  ©)[a£  +  2a^  {coeffw)  +  a^,  +  2^^  {coeffs) 

+  2ap{coefft)  +  ap^p  +®\2afp[{coefw)T^  - {fluxcoefjfw)q'\] 

+  2as  [{coejfs)T^  -  {fluxcoeffs)q'\  ] 

+  2ap  \coefft)T^  -  {fluxcoejft)q" -{radcoefft)Tp  ] 

+  (1  -  ©)[2a^  [(coeffw)T^  -  {fluxcoejfw)q\  ] 

+  2as  (coc#)r„  -  {jluxcoejfs)q'\  ] 

+  2ap  \coejft)T^  -  {fluxcoefft)q'' -(radcoeffi)Tp  ]  (3 .44) 


by  putting  heat  flux  terms  from  equation  (3.7),  equation  (3.42)  may  be  'written  as 

[a“  +  @[ap  +  2a^  {coejfw)  +  +  2a  ^  {coejfs)  +  2ap  {coefft)  +  Og  h 

=  ac  [©T’t  +  0  -  ®)T!  ]+  [®n-  +  (1  -  ®)T°  ]+  a,  [©r,  +  (1  -  e)T‘  ] 

+  [o J  -  (1  -  0)[a J  +  20,  l^coeffic)  +  <i„  +20^  (coeffs) 

+  2ap  {coejft )  +  Jt/  +  \la„  {coeffw)  +20^  {coejfs)  +  2ap  {coejft'^,^ 

-  2ayp  {fluxcoeffw)q'\  -2a  ^  {fluxcoeffs  )q'\  -2a  j  {fluxcoefft  )q" , 
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-  2aj  -2a^  -2aj.  {radcoejf)Tp 


(3.45) 


and,  the  discretised  equation  may  be  written  as 


[a“  +  <d[ap  +  2a^  {coejfw)  ■\-a^  +  20^  {coeffs)  +  20^.  icoejft)+  Og 
-&aj, -&a,T,  =(l-®)[ajr;  +aX  +‘>*7'.“] 

+  a“  -  (1  -  @)[ag  +  2a^  icoeffw)  +  a^  +  2a ^  {coeffs)  +  2^^  {coefft)+  ]]r^° 

+  [2a^  {coejfw)  +  2a ^  {coeffs)+  2aj  {coejft^^  -  2ag,  {fluxcoeffw)q'\ 

-  2a  s  {fluxcoeffs)q'\  -2ap  {jluxcoeffi)q" ,  -2ap  {jluxcoeffi)q'\^^^^^ 

-  2ar  {jluxcoeffi)q'\,umg  {radcoeffjTp  (3.46) 


3.  Front-Left  Edge 


Figure  3.5  Front-Left  Edge 


By  using  equations  (3.28-3.31),  equation  (3.27)  can  be  written  as 


31 


T  -T 

V  ^PE  J 


kA 


T  -T 

KA^ 


t+^\ 


\it+ 


T  -T 

kA^ 

K  ^PN  ) 


T  -T 

kA-^ 

\  ^SP  J 


\it 


/+A/ 
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k,At 


Tr-T, 
& 


PT  J 


k  A 

^b-^b 


\ 


dz 


BP 


\dt 


(3.47) 


In  equation  (3.47),  7”^  and  take  the  value  of  equation  (3.8).  By  using  these  values 
and  equation  (3.35),  we  have 
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(T,-T,) 


+(1-0) 


kM. 


2k Jp  +A^spT^  ^sp 

2k ^ 

f 


+M(r^  -r,) 


PT 


dc 


/J'O  jnO\  ^wA  (  jrO  2^H’7p  Q  >y  3)C^^p 

^  E  pj  c*.  \  P  .  7  0 


PE 


Sc, 


WP  V 


2k yf  '^A^fyp  J 


I  k„A  ,jfi  qrfi^  2k^A^  I  -  2kJI^  q  ^^^p 

■'■  o  v^zv  ^p)  s.  oj  r~c 

^PN  ^sp  \  2k,+K^sp 


k.A, 


kuA^ 


dz 


(r;_r;)-^(r;-r;) 


PT 


dz 


BP 


(3.48) 


Equation  (3.48)  may  be  arranged  as 


r  AF 
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+^[@7;  +(i_e)r;]+^[©r,  +(i-@)r/]+^[©r,  +(i-©)r;] 

2k 


BP 


+ 


+ 


pc^-(l-0) 

At 


& 


2k. 


-  +  • 


^PE  ^wp 


2k ^  +  h^&Cj^  j 


K^r, 

n  n 


PN 


SP  \ 


^K+K^sp) 


+!^+!^\\T^+e\ 

Szpj  Szjgp 


2k^A^ 


Sx 


WP 


h„S)CfypT^  q  ^  Sx^J, 


2k^  +  h^Sx^rp  j 


SP 


2k,A, 


K^spT«  -q”s  ^si-l 

+  (1-0) 

2k 

i 

1 

1 

^WP 

<  2^^  +  h^Sxjyp  J 

SP 


K^spT<.-q'\  ^sp 
2k,+h^^sp  J 


(3.49) 


Eqxiation  (3.49)  may  be  re-arranged  as 


/V*  { 

KA  ,  ^KA^ 

^  h  ^ 

^w^WP 

.  K\  , 

hAs,  1 

T  KJ 
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1  1 

^PN  ^SP 
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^  CP 


AV 
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PE 


2k  A 
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,  KAn  ,  '2-k,A^ 


K^sp 


^PiV  ^SP  '^^s^SP  J 


,  ,  KAj 


&pr  ^BP 


ppO 


+  0 


2k  ^A^ 

h^dx^^T^  q  ^  &Cjyp 

2k,A, 

r  K^spTo. 

SXjYp 

^  2k^  +  h^dxjyp  2k^  +  h„Sx^  j 

"T . 

^SP 

p^s  ^s^SP 

q\^sp  1 

+  (1-0) 

2k  ^A^ 

h^Sxj^T^ 

„ll  ^ 

q  w  ^ivp 

2k^  +  J 

Sx^rp 

[2k„+KSx^ 

2k^+Kdx^j 

^KAs 


SP 


q"s 


K^spT^ 

2k^  2k^+h^dy^  ) 


(3.50) 
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By  putting  the  coefficients  from  the  previous  section,  equation  (3.50)  may  be  written  as 


[flp  +  ©[flf  +  lOfy  (coeffW)  +  2^5  (coeffs)  +  aj  +  ag  Jt,, 

=  +(l-e)rj“  +a„^T„  +(l-0)r,“]+a4®7>  +(1-®)J?] 

+  “a  [®^B  +  (1  “  ®)?» ] +[“’?-(■  -  ®)[‘'b  +  (coefiw)  +o,, 

+  lOs  {coeffs)  +  flr  +  Jr/  +  0[2a^  [{coeffw)T^  -  {jlvxcoejf\^q'\ ] 

+  la^  [{coeffs)T^  -  {fluxcoeffs)q\  J  +  (1  -  ©>[2^;^  [{coeffv)T^  -  {fluxcoeff^q\  ] 

+  los l{coeffs)T„  - {fliacoeffs)q'\  J  (3.5 1) 


Equation  (3.51)  may  be  arranged  as 


+  ©[flg  +  2a^  {coeffw)  +  2a  ^  {coeffs)  +  aj  + a  ^ 

=  a,[®T,  +(i-©)r/]+a^[©r^  +(i-©)r/]+fl,[©rr  +(i-©)r/] 

+  [©r^  +  (1  -  ©)r/  ]+  -  (1  -  ©)[«£  +  2a^  {coeffw)  +0^^ 

+  las  {coeffs)  +  a^.-^  a ^  Jt/  +  \la^  {coeffw)  +  las  {coeffs)]r^ 

-  lafy  {fluxcoeffw)q'\  -las  {fluxcoeffs)q'\  (3.52) 


Finally,  the  discretised  equation  may  be  written  as 


[aj  +  0[oj  +  2a,  (coeffw)  +a„  +  2as  (coeffs)  +  a,  +  a,  JT,  -  ©a^r^  -  @a„r„ 
-0a,r,  -@a,r,  =  (1  -  @)[a,r; + aX + + “.T^h  [4  -0  -0)k 

+  la,y  {coeffw)  +  +  las  {coeffs)  +  aj.  + a ^  Jt/  +  [la^  {coeffw)  +  las  {coeffs)Yoo 
-  la^y  {fluxcoeffw)q'\  -2a ^  {fli4xcoeffs)q'\  (3-53) 


4. 


Bottom  Face 


T 


Figure  3.6  Bottom  Face 


By  using  equation  (3.30-3.33),  equation  (3.29)  may  be  written  as 
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(3.54) 


In  equation  (3.54),  Tg  takes  the  value  of  equation  (3.8).  By  applying  this  and  equation 
(3.35)  to  equation  (3.54),  we  have 
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(3.55) 


Equation  (3.55)  may  be  arranged  as 
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(3.56) 


Equation  (3.56)  may  be  re-arranged  as 


AV 

At 


+  0 


k.A, 

e  e 
Sxpg 


+ 


^  j  ^S-^S 

SXwp  ^PN  ^SP 


k,A, 

+-^-1- 


dz 


PT 


2k, A, 

SZbp 


^  K&bp 

\2k,  +  h,dz^p 


36 


+a-0Fi]+%i=-[®r^  +(i-0)r;]+^[©r„  +a-®)2;“] 

^PE  OXjyp  ^PN 

+-|^[©JV  +o-0)r/]+^[@rr  +(i-©)r;]+U4^-(i-@{%i- 

qy^  oZp^  |_  A?  |_  oXp^ 

I  ^  ^n-^n  I  I  ^t-^t  I  \^BP  jrO  ^  0 

^WP  ^PN  ^SP  ^PT  ^BP  \^^b'^\^Bp)  _  ^BP 

^b^BP^oo _ ^  b  ^BP  ^^b-^b  Bp'^x 

\J2,kb  +hi^&gp  2^J  +hj^&pp  j  ^BP  Kp'^b  ^b^BP 


<l"b  ^BP 

2ki,  +\dzpp  j_ 


(3.57) 


By  using  the  defined  coefficients,  equation  (3.57)  may  be  written  as 


[a°  +  @[ap  +  +aj^  +  +  Oj.  +  2a  p  {coeffl})^p  =  [@Tp  +  (1  -  0)7’^°  ] 

+a^[0r^  +(i-0)r^]+a^[0r^  +(i-0)r^°]+fl,[0r,  +(i-0)r/] 

+  Qp [0rr  +  (1  - &)Tp  ]+  [a^  - (1  - 0)[a£  +  %  +as+ap+  2a p {coeffl))^p 

+  ®]^ap  [((coe^)r„  )  -  {{fluxcoefflj)q'\  )] 

+  (1  -  0)[2aB  [((coc^)r„  )  -  {{flvxcoeffli)q'\  )]  (3.58) 


And,  the  discretised  equation  may  be  written  as 


[al+e[ap  +a^  +%  +^5  +ap  +  2ap{coeffb)^p  -©a^Tp  -©a^pTyp  -©a^^T^ 

-  ©apTp  -  ©apTp  =  (1  -  0)[«£r°  +  ajpT°  +  aj^T^  +  apT^  +  a^Tj  ]+  [a® 

-  (1  -  ©)[ap  ■¥a^+a^+as+ap  +  2a  p  {coeffi>)^p  +  \2ap  {coeffl})]r^ 

-  2a  p  {fltacoeffli)q'\  (3.59) 


The  resulting  discretised  equations  for  the  other  parts  of  the  workpiece  are  as  follows 
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5.  Left-Back-Top  Corner 


\a°  +  ©[«£  +  2a  ff,  (coe^fw)  +  2a  {coeffri)  +  a^  +  2a^  {coefft)  +  Jr,  -(da^T^ 

- =  (1  - +  cisTs  +  ]+  [ap  -  (1  -  ®)["£  +  2a^ (coefjw) 

+  2a  {coeffii)  +  a^  +  2aj.  {coefft)  +  a^  ]r°  +  [2a„-  (coeffw)  +2af^,  (coeffn) 

+  2a j. {coefft)]r^  - 2a„, {fluxcoeffw)q'\ -2aj^{fluxceeffri)q'\ -2a j {fluxcoefft)q'\ 

-  2aj-  {fluxcoefff)q'\^^^  -2a  j  {fluxcoefft)q\,,,„^ 

-  2a  J  {radcoeff)Tp  (3 .60) 


6.  Right-Back-Top  Comer 

[aj  +  ©[2a£  {coeffe)  +  a^  +  2a f,  {coeffh)  +  a^  +  2a^  {coefft)  +  -  ©a^T^ 

-@affs -®affp={\-Q)[a^T^  +a5r;  +affl]+[al -{\-@)[2ap{coeffe) 

+  %  +  '2-0  s  {coeffh)  +  as+  2aj.  {coefft)  +  a^  +  [20^  {coeffe)  +2a^  {coeffh) 

+  2ap  {coefft)]r„  -  2a ^  {fluxcoeffe)q'\  -2a ^  {fluxcoeffn)q'\  -2a j.  {fluxcoefft)q'\ 

-  2ap  {fluxcoefft)q'\^^^^  -2ap  {fluxcoefft)q 

boiling 

-  2aj.  {radcoeff)Tp  (3.61) 


7.  Right-Front-Top  Comer 

[a;+©[2af  {coeffe)  +  a^  +  a^  +  2a  ^  {coeffs)  +  2a  j  {coefft)  +  a^  '^p-QayyT^ 

-  ©afff,  -  @agTp  =  (1  - ®)[a^T°  +  a^T^  +  a^Tp  ]+  -  (1  -  ©)[2a£  {coeffe) 

+  a^  +  a^  +  2a p  {coeffs)  +  2ap  {coefft)  +  a^  +  [2a^  {coeffe)  +2ap  {coeffs) 

2ap  {coeffi)]r^  -  2a ^  {fluxcoeffe)q'\  -2a ^  {fluxcoeffs)q'\  -2aj  {fluxcoefft)q'\ 

-  2aj  {fluxcoeffi)q'\^^„-2ap  {fluxcoeffi)q 

boiling 

-  2aj{radcoeff)Tp 


(3.62) 


8. 


Left-Back-Bottom  Corner 


[aj  +  ©[«£  +  2ajy  {coejfw)  +  2a^{coeffit)  +  as+aj+  2aB(coeffl))^p  - 
- - &aj.Tj  =  (1  - &)[apTp  +  0^1°  +  aj.Tj\+  [a°  - (1  - 0)[o^  +  2a^ {coeffw) 
■\-2aj^  (coeffn)  +  as+aj+  2apicoeffb)Yp 

+  [2a^ (coeffw)  +2af^(coeffn)  +  2ap {coeffb)^„  - 2a^ (flt4Xcoeffw)q'\ 

-2aj^  {flvxcoeffn  )q'\  -2a p  (^flvxcoeffb  )q'\  (3.63) 


9.  Right-Back-Bottom  Comer 

[otj  +  0[2a^  {coeffe)  +  ayp+  la^  {coeffn)  +  ^5+0^.+  2a ^  {coeffb)^p  -  ®a^T^ 

- @asTs  - ®apTj  =  (1  - ©)[a^r^  +  aJs  +  «r^r  ]+  Vl  “ 0  “ ©)[2a£(coe#e) 

+  afp+2apf  {coeffn)  +  +  a^.  +  2a  p  {coeffb^l 

+  \lap {coeffe)  +20^^ {coeffh)+  2aB{coeffb)Y'^  -  2a ^ {fluxcoeffe)q'\ 

-  2a  ^  {ff.ux:coeffh)q'\  -2a  p  {fluxcoeffb)q'\  (3.64) 


10.  Left-Front-Bottom  Comer 

[aj  +  @[ap  +  2a„r  {coeffw)  +  aj^+2as  {coeffs)  +  +  2a p  {coeffb)^p  -  ®apTp 

-  ®^nTn  -  =  (1  -  ®)[a£7’£  +  +  ^tTt  ]+[«“-  (1  -  ®)[«£  +  '^w  {coeffw) 

+  +  2as  {coeffs)  +  a^.  +  2a  p  {coeffb)^p 

+  [2a^  {coeffw)  +2ap  {coeffs)  +  2ap  {coeffb)^„  -  2a^  {fl.tacoeffw)q'\ 

-2as  {fluxcoeffs)q'\  -2a  p  {flvxcoeffb)q'\  (3 .65) 


11.  Right-Front-Bottom  Comer 

[al+®\2ap  {coeffe)  +  affr  +  Off  +  2a ^ {coeffs)  +  +  2ap{coeffb^p  - ®a^T^ 
-@a^T^  =(l-0)[a^7;“  +a^T° +apTp^]+[al -{l-@)l2ap{coeffe) 

+  ajp  +  a  j^+  2a  p  {coeffs)  +  +  2a  p  {coeffb)^p 

+  \2ap {coeffe)  +2^^ {coeffs)  +  2ap {coeffb)^„  -2a p {flvxcoeffe)q'\ 

-  2as  {flwccoeffs)q'\  +2ap  {fluxcoeffb)q'\  (3.66) 
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12.  Top  Face 


-@aj, -@aj, -®a,T,={\-®)\aJI+aX'+ciX  +  cisT^s+ci,Tl] 

+  Wp  “  (1  -  ®)[^£  +  %  +  ^N  +  ^5  +  2ar  (.coeffi)  +  a^  Jt/ 

+  \laj{coefft)]r^  -  2a j. {fluxcoeffi)q'\ -la^. {fluxcoefft)q'\^^^ 

-  laj.  (fluxcoeffi)q”^„,„^  -la^-  {radcoejfX  (3 .67) 


13.  Front  Face 

[op  +  ©[flf  +aiy  +  a  +  2a^  {coeffs)  +  a.j-+a^  -&aX 
-  -  &a,T,  -  ea,T,  =  (I  -  ©)[<7jrj"  +  a,.T°  +  aX  +  “X  +  o.!",”] 

+  [aj-(l-©)[a^+a^  +  +  2a  ^  {coeffs)  +  a  j-  +ag  fc' 

+  [2a s  (coeffs)y„  -  2a ^  {fluxcoeffs)q'\  (3 .68) 


14.  Back  Face 

[lOp  +  ©[^f  +  Off,  +  2a {coeffn)  +  a^  +  a-j-  +  a^  ~^ffp  ~  Qaff^  —  ©dt^T^ 

-©a^Ts  -©off, -©aj,  =(l-©)[a^7’;  +aX +^3^° +aX +(‘bTb] 

+  [a^  -  (1  -  ©)[fl£  +  aff,  +  2a  ff  {coeffn)  +  Ih" 

+  [2a^  {coeffh)y^  -  2a ^  {fluxcoeffh)q'\  (3 .69) 


15.  Left-Lateral-Face 

Ifip  +  ©[«£  +  2aff,  {coeffw)  +  -  ©affp  -  ©afff^, 

-&asT, -earTr  -©a^r,  =(l-®)[ajr/+<i,r;  +aX*oX  ^a.^S] 

+  [a"-(l-©)[<ij+2a,  {coeffw)  +  +  fltj.  +  Ih" 

+  [2%  {coeffw)]r^  -  2afy  {fluxcoeffw)q'\  (3.70) 
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16.  Right-Lateral-Face 


[al  +  0[2fl^ (coeffe)  +  -  @a^T^  -  &aj^Tj^ 

—  ®a^Tg  —©Oj-Tj-  —&cigTg  =  (1  — 0)[a(^7’^  +aj^Tj^  +cigTg  +cij-Tj-  +^^£75] 

+[ii;-(i-0)[2<jj  (coeffe)  +  <7^  +  (3^  +  +  flj-  +  <3^  fc 

+  \laE{coejfe)%  - 2a ^ {fluxcoeffe)q'\  (3.71) 


17.  Front-Top  Edge 

[al  +  0[a£  +afy  +  a j^+  la^  {coeffs)  +  la^  (coefft)  +  -  ©a^TE 

—  &ajfrTjy  —  ©a^fTfi  —  ©a^Tg  =  (1  ~  +  ^w'^w  ^n'^n  ] 

+  [al  -  (1  -  0)[fl[jj  +afy+a^+  2a ^  (coeffs)  +  2a j.  (coeffi)  +  ijr^ 

+  \2as  (coeffs)  +  20^  (coeffi)]r„  -  2a g  (Jlvxcoeffs)q'\  -2a j.  (fluxcoeffi)q'\ 

-  2a J  (fluxcoeffi)q'\,,,^  -2aj.  (fluxcoefft)q” 

-2aj(radcoeff)T^  (3.72) 


18.  Front-Bottom  Edge 

[al  +  @[ag  +  <7^  +  +  2as (coeffs)  +  aj+  2ag(coeffb)^p  - ea^Tg 

—  ®ajfT If  —  —  0<3j.7^  ~  (1  ~  ^w'^w  ^n'^n  ] 

+  [al  -  (1  -  0)[a£  +  a^  +  a f^+  2a g  (coeffs)  +  Oj.  +  2a g  (coeffb)^^ 

+  [2as  (coeffs)  +  2a  g  (coeffb)]r^ 

-  2a  g  (fluxcoeffs)q'\  -2a  g  (fluxcoeffb)q'\  (3 .73) 


19.  Back-Top  Edge 

[al  +  @[ag  +ajy+  2a  f,  (coeffh)  +  +  2ap  (coeffi)  +  a  g'^p-  QagTg 

-®a^Tjp  -©aglg  -@agTg  =(l-0)[a£7’£  +a^T°  +agTs  +agTg] 

+  [al (1  -  0)[fl^  +afp  +  2a (coeffh)  +  +  2ap  (coeffi)  +  Og  Jt/ 

+  \2af, (coeffh)  +  2ap (coeffi'^^  - 2a^(fluxcoeffh)q"„ -2% (flwccoeffi)q" , 

-  2ap  (fluxcoefft)q”,^„,-2ap(fli4XCoeffi)q"g^ai„g  -2ap(radcoeff)Tp  (3.74) 
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20.  Back-Bottom  Edge 


[a®  +  ©[«£  2a {coeffii)  +  0^+0^+  2a ^  {coeffl})^p  -  Qa^T^ 

-ea^T^ -@aj, -®apTp  =(l-e)[apT°  +  a^T^  +a,r;  +fl,r;] 
+  [aj  -  (1  -  ©)[^^£  +  +  2oa-  (coeffn)  +  0^+0^+  2a ^  (coeffly)^p 

+  l2ai^  {coeffn)  +  2a  p  {coeffb)]r^ 

-  2a  ^  {fluxcoeffn)q'\  -2a  p  {fluxcoeffb)q'\ 


21.  Front-Right  Edge 

[flf®  +  ®\2ap  {coeffe)  +  +  205  {coeffs)  +  ih  “ 

-@aj^  -®apTp  -QapTp  =(l-©)[a^r^  +a^T^^  +apT^  +  a pT^] 
+  [«J  -  (1  -  ©)[2a£  (coe#e)  +  +  2^^  {coeffs)  + 

+  [2a£  {coeffe)  +  2a  ^  {coeffs)]r^ 

-  2ap{fliixcoeffe)q'\  -2a  ^  {fluxcoeffs)q'\ 


22.  Back-Left  Edge 

[a J  +  ©[fl^  +  2ajp {coeffw)  +  2a ^ {coeffn)  +  a^  +  a^  +  Op'^p  - ©OpTp 

-  eapTp  -  QapTp  -  BopTp  =  (1  -  ©)[a^r;  +  a.r/  +  a^T^  +  OpT^  ] 

+  [ofj  -  (1  -  ©)[«£  +  2apr  {coeffw)  +  2a ^  {coeffn)  +  0^+0^+^^  ik* 

+  [2a^  {coeffw)  +  2a  j^{coeffh)]r^ 

-  2a^  {fluxcoeffw)q'\  -2a  {fluxcoeffh)q'\ 


23.  Back-Right  Edge 

[aj+©[2a^  {coeffe)  +  +  2a  ff  {coeffn)  +  +  atj-  +  -©%r^ 

-  ®asTs  -  ®apTp  -  QopTp  =  (1  -  ®)[a^T^  +  a^T^  +  a^T^  +  OpTp  ] 

+  [aj  -  (1  -  0)[2a£  {coeffe)  +  a^+  2a ^  {coeffn)  +  +  Oj.  +  Jt^® 

+  ]?Op  {coeffe)  +  2a  {coeffn)]r^ 

-  2a  p  {fluxcoeffe)q'\  -2a  {fluxcoeffh)q'\ 


(3.75) 


(3.76) 


(3.77) 


(3.78) 
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24.  Left-Lateral-Top  Edge 

\a°p  +  ©[«£  +  lojy  (coeffw)  +  +  la^  (coefft)  +  ag  fc  -  &apTp 

-  @a,T,  -  ©aJs  -  =  (1  -  @)[apT°  +  aX  +  +  ^bTb  ] 

+  [fl°  - (1  - ©)[ag  +  2a^ (coeffw)  +  af^+as+2ap (coefft)  +  jr^ 

+  [2a^  (coeffw)  +  la^.  (coefft)]r„  -  (fluxcoeffw)q'\  -la^  (jluxcoefft)q'\ 

-  lap  (flvxcoeffi)q'\^^-2ap  (fluxcoeffi)q" 

boiling 

-lap  (radcoeff)Tp  (3.79) 


25.  Left-Lateral-Bottom  Edge 

[aj  +  ©[a^  +  lajy  (coeffw)  +  +  2a5  (coeffb)^p  -  ©a^Tg 

-  @a,T^  -  ©a,r,  -  (dapTp  =  (1  -  ©)[a^r;  +  aX  +  ^X  +  cipT^  ] 

+  -  (1  -  ©)[«£  +  2a^  (coeffw)  +  +  2ag  (coeffb)^p 

+  \la^  (coeffw)  +  lap  (coeffb)]r„ 

-  lOfp  (fluxcoeffw)q'\  -2aB(fluxcoeffb)q'\  (3.80) 


26.  Right-Lateral-Top  Edge 

[a°  +  ©[2a£  (coeffe)  +  +  lap  (coefft)  +  a^  Ih  -  &ajfrTjp 

-  ®aX  -  ®^sTs  -  ®cipTp  =  (1  -  ®)[aX  +  ciX  +  ^sTs  +  ] 

+  [a  J  -  (1  -  ©)[2a£  (coeffe)  +  +  lap  (coefft)  +  a^  Jt/ 

+  \lap  (coeffe)  +  lap  (coeffi)]r^  -  la^  (fluxcoeffe)q'\  -lap  (flwccoeffi)q'\ 

-  lap  (fluxcoeffi)q'\^^-lap  (flvxcoeffi)q 

boiling 

-lap(radcoeff)Tp  (3-81) 
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27.  Right-Lateral-Bottom  Edge 


[al  +  @\la^  {coejfe)  +  +  a,.  +  2ag{coeffl))^- 

-  -  ®asTs  -  ©arT,  =  (1  -  @)[a^K  +  ciX  +  ] 

+  [aj  -  (1  -  ©)[2a^  (coe#e)  +  +  a^.  +^5+0;.+  2a ^  {coeffl})^^ 

+  [2a£  (coe#e)  +  2a  ^  {coeffl})]r^ 

-  2a  E  {fluxcoeffe)q'\  -2a  ^  {jluxcoeffb)q'\ 


(3.82) 
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IV.  RESULTS  AND  DISCUSSION 


Results  are  presented  for  the  different  cases  discussed  in  greater  detail  below.  As 
noted  earlier,  a  variable  sized  and  moving  numerical  mesh  was  used  in  such  a  way  that 
the  arc  was  always  positioned  at  the  center  of  the  mesh  where  the  spacing  is  the  finest. 
The  goal  is  to  be  able  to  resolve  the  large  temperature  gradient  features  around  the  arc 
and  yet  not  incur  a  large  overhead  of  computer  resources,  which  would  be  required,  if  the 
grid  was  uniformly  fine  all  over.  Tbis  strategy  however  did  require  that  a  separate  mesh 
generating  rotrtine  be  used  which  did  demand  some  extra  computational  resources.  The 
weld  pool  region  was  also  modeled  as  a  solid  region  but  with  a  thermal  conductivity 
higher  than  the  surrounding  unmelted  region  to  simulate  the  effects  of  weld  pool 
convection.  The  discontinuity  in  the  thermal  conductivity  boundaries  was  handled  using 
the  standard  technique  of  employing  harmonic  averaging  at  the  boundary.  Since  the 
coefficients  of  the  system  of  equations  depend  on  the  temperature,  an  iterative  solution 
technique  was  used  to  achieve  convergence  in  such  a  way  that  the  maximum  temperature 
difference  between  two  consecutive  iterations  at  any  grid  point  was  no  more  than  0.1  °C . 

The  numerical  solution  method  was  used  to  examine  different  cases  in  fireshwater 
for  a  40-mm-thick  70  x  90  mm  worlq)iece  with  a  moving  heat  source  in  the  positive  y- 
direction. 

Case  la: 

The  validity  of  the  numerical  model  was  compared  to  Rosenthal’s  three- 
dimensional  solution  for  a  moving  heat  source.  At  this  point,  convective,  radiative  and 
boiling  surface  thermal  conditions  were  not  considered.  A  constant  thermal  conductivity 
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value  and  a  point  heat  source  were  used  in  the  calculations.  The  input  data  used  in 
computations  for  case  la  are  shown  in  Table  1.  The  calculations  were  made  up  to  3.6 
seconds.  The  numerical  results  show  excellent  agreement  with  the  analytical  results  of 
Rosenthal  (Figure  4.1  and  Figure  4.2). 


Workpiece  :  Length  =  70  mm,  Width  =  90mm,  Thickness  =  40mm 
Water  temperature:  T  =  27°C  (300.15  K) 

Power  input  into  workpiece  :  Q  =  2544  W 
Arc  torch  speed  :  Vy=  4  mm/s 
Radius  of  heat  input  distribution  :  ro  =  4.5  mm 
Thermal  conductivity  :  k  =  53  W/m  K 
Density  of  steel :  p  =  7854  kg/rr^ 

Specific  heat  of  steel :  Cp  =  509.3  J/kgK 


Table  1.  The  input  data  used  in  computations  for  case  la 
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(b) 

Figure  4.2  The  numerical  model  point  heat  source  solution: 

(a)  Profile;  (b)  Front-view. 

Case  lb  : 

In  this  case,  the  surface  thermal  boundary  conditions  of  boiling  heat  transfer  were 
now  included.  Other  relevant  data  is  given  in  Table  2. 

Workpiece  :  Length  =  70  mm,  Width  =  90mm,  Thickness  =  40mm 
Water  temperature:  T  =  27‘’C  (300.15  K) 

Saturation  temperature:  T=  100  °C  (373. 15  K) 

Water  depth  :  1  =  Oft 

Total  pressure:  P  =  101.325  kPa 

Power  input  into  workpiece  :  Q  =  2544  W 

Arc  torch  speed  :  Vy  =  4  mm/s 

Radius  of  heat  input  distribution  :  To  =  4.5  mm 

Thermal  conductivity ;  k  (W/m  K) 

53  -  0.04  (T-300)  300  <  T(K)  <  1000 

25  +  6.25  X  ia^(T-1000)  1000  <  T(K)  <  1800 

125  T(K)  >  1800 

Emissivity  :  S  =  0.82 _ 

Table  2.  The  input  data  used  in  computations  for  case  lb 
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The  top  surface  temperature  values  at  0.5  seconds  can  be  seen  in  Figure  4.3  and 
Figure  4.4.  At  0.5  seconds,  the  temperature  distribution  is  almost  symmetric  about  the 
position  of  the  arc.  The  temperature  distribution  around  the  arc  center  is  above  the 
melting  temperature  of  the  worlq)iece  (Tm  =  1800  K).  To  see  the  depth  of  the  weld  pool 
penetration,  the  melting  temperature  contoms  were  plotted  for  different  surfaces  (Figure 
4.5)  and  the  weld  pool  depth  was  shown  through  the  melting  temperature  points  by  using 
a  curve-fit  (Figure  4.6).  Because  of  the  brief  reaction  time,  the  arc  heat  input  penetration 
to  the  worlqpiece  is  limited  and  a  well-formed  weld  pool  cannot  be  seen.  [Ref.  16] 


Figure  4.3  Top  surface  t=0.5  sec 
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3000 


y-axis  (m) 


Figure  4.4  Top  surface  (profile)  t=0.5  sec 


Figure  4.5  The  weld  pool  surface  characteristics 
by  the  T=1800  K  contours  t=0.5  sec 
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Figure  4.6  The  weld  pool  (front-view)  t=0.5  sec 


At  t  =  2.0  sec.,  the  results  reveal  that  the  temperature  distribution  has  already 
reached  steady  state.  Due  to  the  moving  heat  source,  the  small  change  of  slope  of  the 
temperature  distribution  curve  can  be  observed  clearly  (Figure  4.7  and  Figure  4.8). 
Because  of  the  increased  reaction  time,  the  melting  temperature  contours  are  seen  till  the 
fourth  surface  from  the  top  (Figure  4.9).  The  weld  pool  width  increased  to  8  mm  and  its 
depth  to  2.1  mm  (Figure  4.10). 


51 


Figure  4.7  Top  surface  t=2.0sec 


Figure  4.8  Top  surface  (profile)  t=2.0  sec 
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Figure  4.9  The  weld  pool  surface  characteristics 
by  the  T=1800  K  contours  t=2.0  sec 
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Figure  4.10  The  weld  pool  (front-view)  t=2.0  sec 
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At  t=  4.5  sec.,  the  change  of  slope  of  the  temperature  distribution  behind  the  arc  is 
more  evident  in  Figure  4.1 1  and  Figure  4.12.  The  melting  temperature  contours  are  still 
seen  till  the  fourth  surface  from  the  top  (Figure  4.13).  But,  the  weld  pool  depth  reached  to 
3  mm  with  a  width  of  8  mm  (Figure  4.14).  These  results  show  a  good  agreement  with  the 
data  from  Ref.  16. 


Figure  4.11  Top  surface  t=4.5  sec 
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Figure  4.13  The  weld  pool  surface  characteristics 
by  the  T=1800  K  contours  t=4.5  sec 
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Figure  4.14  The  weld  pool  (front-view)  t=4.5sec 


Case  2, 3, 4  : 

In  the  following  cases  we  do  not  consider  the  temperature  distribution.  Instead, 

the  cooling  time  that  elapses  between  800  °C  and  500  °C  is  examined  for  freshwater  at 
different  temperatures  and  welding  depths.  The  input  data  used  in  calculations  are  shown 
in  Table  3,  Table  4  and  Table  5. 
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Workpiece  :  Length  =  70  mm,  Width  =  90mm,  Thickness  =  40mm 
Water  temperature;  T  =  3  °C  (276. 15  K) 

Saturation  temperature  ;  T  =  114.88  °C  (388.03  K) 

Water  depth  :  1  =  22  ft  (6. 706  m) 

Total  Pressure  :  P  =  168. 75  kPa 

Power  input  into  workpiece  :  Q  =  3900  W 

Arc  torch  speed :  Vy  =  2. 75  mm/s 

Radius  of  heat  input  distribution  :  ro  =  4.5  mm 

Thermal  conductivity ;  k(W/mK) 

53  -  0. 04  (T-300)  300  <  T(K)  <  1000 

25  +  6.25  X  m^CT-lOOO)  1000  <  T(K)  <  1800 
125  T(K)  >  1800 

Emissivity :  S  =  0.82 _ 


Table  3.  The  input  data  of  calculations  at  T  =  3  C  and  1  =  22  ft 


Workpiece  :  Length  =  70  mm.  Width  =  90mm,  Thickness  =  40mm 
Water  temperature:  T=10°C  (283.15  K) 

Saturation  temperature  :  T=  112.62  °C  (385.77 K) 

Water  depth  :  I  =  18  ft  (5.486  m) 

Total  Pressure  :  P  =  156.492  kPa 

Power  input  into  workpiece  ;  Q  =  3900  W 

Arc  torch  speed :  Vy  =  2. 75  mm/s 

Radius  of  heat  input  distribution  :  ro  =  4.5  mm 

Thermal  conductivity ;  k  (W/m  K) 

53  -  0. 04  (T-300)  300  <  T(K)  <  1000 

25  +  6.25  X  10-^(7-1000)  1000  <  T(K)  <  1800 

125  T(K)  >  1800 

Emissivity :  6  =  0.82 _ 


Table  4.  The  input  data  of  calculations  at  T  =  10  C  and  1  =  18  ft. 
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Worlq)iece  :  Length  =  70  mm,  Width  -  90mm,  Thickness  =  40mm 
Water  temperature:  T  =  31°C  (304.15  K) 

Saturation  temperature  :  T  =  115.44  °C  (388.59  K) 

Water  depth  :  1  =  24  ft  (7.315  m) 

Total  Pressure  :  P  =  171. 762  kPa 

Power  input  into  workpiece  :  Q  =  3900  W 

Arc  torch  speed  :  Vy  =  2. 75  mm/s 

Radius  of  heat  input  distribution  :  ro  =  4.5  mm 

Thermal  conductivity ;  k  (W/m  K) 

53  -  0.04  (T-300)  300  <  T(K)  <  1000 

25  +  6.25  X  10-^(7-1000)  1000  <  T(K)  <  1800 

125  T(K)  >  1800 

Emissivity :  £  =  0.82 _ 


Table  5.  The  input  data  of  calculations  at  T  =  31  C  and  1  =  24  ft. 


The  resulting  graphs  of  the  water  temperature  and  welding  depth  effects  on  the 

peak  temperature  and  the  cooling  time  (between  800  °C  and  500  "C )  for  a  point  initially 
1.25  mm.  behind  the  arc  heat  soiuce  are  shown  in  Figure  4.15,  Figure  4.16  and  Figure 

4.17.  The  cooling  times  for  800  -  500  °C  temperature  range  are  0.19  seconds(Twater=3 

‘’C,  1=22  ft.),  0.20  seconds  (Twater=10  °C ,  1=18  ft.)  and  0.33  seconds  (Twater=31  °C ,  1=24 
ft.)  respectively.  These  cooling  times  are  significantly  different  firom  those  in  the  previous 
studies  in  the  literature.  Based  on  the  models  in  this  case,  it  is  seen  that  boiling 

phenomena  must  be  considered  for  surface  heat  losses  from  the  weld  metal.  For  800  °C 

to  500  °C  range  of  the  weld  metal,  film  boiling  exists  and  the  surface  is  completely 
covered  by  a  vapor  blanket.  A  high  amount  of  heat  transfer  firom  the  weld  metal  to  the 
water  environment  occurs  by  conduction  across  the  vapor  film.  The  other  reason  is  that 

due  to  the  thermal  conductivity  equations  used  in  the  model,  for  800  °C  to  500  °C 
ranges,  thermal  conductivity  of  steel  is  low.  This  causes  a  decrease  in  conduction  heat 
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transfer  from  the  arc  to  the  studied  point.  But  the  hi^  convection  heat  loss  at  this  point 
causes  a  rapid  cooling  rate.  Another  possible  reason  is  that  because  of  the  low  welding 
speed  and  the  high  thermal  conductivity  in  the  weld  pool,  the  workpiece,  which  has  a 
large  diermal  capacity,  absorbs  most  of  the  heat  energy.  After  the  arc  has  passed, 
temperature  and  thermal  conductivity  at  that  point  drops.  Due  to  the  low  thermal 
conductivity,  the  transfer  of  the  absorbed  heat  energy  from  the  inner  part  of  the 
worlq)iece  to  the  surface  is  not  sufficient  compared  to  the  surface  heat  loss.  This  causes  a 
rapid  cooling  rate  on  the  surface. 

During  welding,  the  high  heat  input  to  the  workpiece  increases  the  degree  of  grain 
growth.  The  grain  coarsening  effect  causes  a  decrease  in  the  grain  boundary  area. 
Because  of  the  large  grain  size  and  the  very  high  cooling  rates,  the  resulting  phase  is  hard 
with  Vi  >  450,  but  with  a  brittle  martensitic  structure. 


Figure  4.15  Cooling  time  and  peak  temperature  for  a  point 
initially  1.25  mm  behind  the  arc  heat  source. 
(Twater=  3  C,  depth  =  22  ft.) 
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Figure  4.16  Cooling  time  and  peak  temperature  for  a  point 
initially  1.25  mm  behind  the  arc  heat  source. 
(Twater  ~  10  Cj  depth  —  18  ft.) 


Figure  4.17  Cooling  time  and  peak  temperature  for  a  point 
initially  1.25  mm  behind  the  arc  heat  source. 
(Twater=  31  C,  depth  =  24  ft.) 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


Compared  to  experimental  methods,  numerical  solution  techniques  are  very  fast, 
effective  and  economical  alternatives.  Therefore,  the  aim  of  this  study  was  to  develop  a 
general  numerical  model  for  transient,  three-dimensional  conduction  heat  transfer 
phenomena  in  underwater  welding  process  on  a  thick  rectangular  plate.  Computations 
were  presented  for  different  cases.  Comparisons  of  current  predictions  with  results  in  the 
literature  showed  good  agreement  and  validated  the  model.  If  the  material,  environment 
and  the  arc  source  properties  are  known,  this  program  can  be  applied  to  the  different 
types  of  metals  under  the  wet  welding  process. 

The  numerical  model  gave  important  temperature-time  data  in  the  critical  HAZ, 
which  in  turn  determines  material  structure.  The  model  also  helped  to  make  prediction  of 
size  of  weld  pool  and  its  evolution  with  time. 

The  principal  recommendations  for  the  future  studies  may  be  summarized  as 
follows. 

1.  The  weld  pool  region  itself  was  modeled,  as  a  solid  region.  But,  owing  to 
the  presence  of  the  temperature  values  higher  than  the  melting  temperature 
in  the  weld  pool,  the  model  must  be  considered  as  a  liquid  weld  pool  with 
the  surroimding  unmelted  solid  region.  So,  the  numerical  model  needs  to 
be  modified  to  account  for  weld  pool  convection  with  the  buoyancy  force, 
the  electromagnetic  force  and  the  surface  tension  gradient  at  the  weld  pool 
surface. 

2.  The  chemical  effects  due  to  interaction  of  the  arc  with  the  surrounding 
water  environment  may  need  to  be  included. 

3.  In  the  model,  boiling  heat  transfer  phenomena  was  accormted  for  a  solid 
surface.  For  a  molten  pool  region,  the  role  of  boiling  on  a  liquid  substrate 
is  unknown  and  needs  to  be  solved. 
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APPENDIX  A.  PROGRAM  STRUCTURE 


In  the  main  program  <aathl0.m>,  non-uniform  grid  spacing  was  constructed  by 
using  the  exponential  function  y  =  he~‘’^  where  b  is  the  coefficient  which  affects  the  grid 
spacing  distribution  and  h  is  the  sum  of  the  grid  spacings  when  x  ^  oo .  The  application 
of  grid  spacing  to  the  workpiece  was  shown  in  Figure  A1 . 


(b) 

Figure  A1  Applying  the  grid  spacing  to  the  workpiece  by  using 
the  exponential  function 
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The  distance  in  front  of  the  moving  heat  source  was  called  as  “front”  and  the 
distance  behind  the  heat  source  was  called  as  “back”.  The  width  and  the  thickness  of  the 
plate  were  defined  with  the  same  terms.  The  variable  “number”  was  used  to  define  the 
number  of  grid  points  in  the  “back”  region.  Variables  “b”  and  “number”  can  be  changed 
to  control  the  resolution  of  the  grid  spacing. 

To  simulate  the  heat  input  from  the  arc  to  the  top  svuface,  it  was  assumed  that  the 
arc  source  heat  flux  distribution  had  a  Gaussian  distribution  on  the  top  surface.  The 
applied  arc  heat  source  was  defined  as  a  matrix,  which  its  size  was  equal  to  the  size  of  the 
top  surface.  The  value  of  heat  flux  at  each  grid  point  was  calculated  by  using  the  x- 
coordinate  and  the  y-coordinate  of  that  point.  The  resulting  arc  heat  source  matrix  was 
applied  to  the  top  surface. 

The  discretised  equations  was  represented  by  the  matrix  equation  [A]  X  =b 
where,  [A]  is  the  coefficient  matrix,  X  is  the  column  matrix  of  the  unknown  temperature 
values  of  the  grid  points  and  b  is  the  column  matrix  of  the  constants.  In  the  coefficient 
matrix,  the  coefficients  of  the  studied  grid  and  the  neighboring  grids  were  written  to  the 
same  row.  An  example  for  5x5x2  workpiece  was  shown  in  Figure  A2.  In  this  example, 
the  coefficient  matrix  may  be  written  as 

cl  c2  0  0  0  c6  0...  c26... 

cl  c2  c3  0  0  0  c7  0...  c27... 

0  c2  c3  c4  0  0  0  0...  c28... 

0  0...  c25...  c45...  c49  c50 
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Figure  A2  The  example  workpiece 


Even  the  row  number,  the  column  number  and  the  surface  number  were  increased, 
the  maximum  number  of  coefficients  for  each  row  does  not  exceed  7.  This  is  for  an 
interior  grid  point  of  a  workpiece  with  a  surface  number  3  or  more.  But,  the  increasing 
size  of  the  coefficient  matrix  increases  the  computing  time  and  the  required  computing 
memory  and  storage  capacity.  To  prevent  this  problem,  tridiagonal  part  of  the  general 
matrix  (Figure  A3)  which  contains  Tw  ,Tp  and  Tp  coefficients  was  formed  as  a  new  ax3 

matrix  (a  =  row  x  column  x  surface).  The  coefficient  values  and  TV.  r^,  Tp  and  Tp 

were  also  formed  separately  and  called  as  <north>,  <south>,  <top>  and  <bottom> 
matrices  respectively.  The  size  of  each  matrix  is  a  x  1  (a  =  row  x  column  x  surface).  The 
constant  values  in  the  right  side  of  the  equation  were  called  as  <right>  matrix. 
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Figure  A3  The  coefficient  matrix 


Beginning  from  the  first  surface,  the  fimction  program  <aasolvelO>  solves  the 
matrix  produced  by  the  main  program  by  using  row-by-row  iterative  sweeping  method. 
The  representation  of  row-by-row  method  was  shown  in  Figure  A4  [Ref.  21,22].  The 
statements  <indexl>  and  <index2>  defines  the  first  and  the  last  element  of  each  row. 


•  Points  at  which  values  are  calculated 

o  Points  at  which  values  are  considered 
to  be  temporarily  known 

■  Known  boundary  values 
— ►  Sweeping  direction 

Figure  A4  Row-by-row  method  representation 


The  matrix  <temp_old>  represents  the  old  temperature  values  of  the  grid  points. 
At  first,  it  was  initialized  and  then  it  took  the  values  of  the  previous  iteration.  The  matrix 
<temp>  represents  the  temperature  values  of  the  grid  points,  which  were  found  from  the 
present  iteration.  The  sizes  of  <temp_old>  and  <temp>  are  (a  x  b)  x  c,  where  a,  b  and  c 
are  the  total  row  number,  colimm  number  and  surface  number  respectively.  The  statement 
<delta_iteration>  represents  the  absolute  value  of  the  difference  between  <temp_old>  and 
<temp>.  If  <delta_iteration>  is  less  0.1,  the  resulting  temperature  value  is  satisfactory. 
The  statement  <count_iteration>  counts  the  iteration  number.  If  the  program  can  not 
converge  after  10  iterations,  <delta_iteration>  increases  0.1  and  continues  to  increase  0.1 
at  every  following  5  iterations  until  the  program  find  a  converged  temperature  value.  The 
radiation  heat  flux  matrix  <radiationl>  and  the  boiling  heat  flux  matrix  <q_boilmg_fi'ee> 
are  only  used  to  calculate  the  heat  flux  from  the  top  surface. .  The  sizes  of  <radiationl> 
and  <q_boilmg_fi:ee>  are  (a  x  b)  x  1,  where  a,  b  are  the  total  row  number  and  col;imn 
number  respectively. 

In  the  <aasolvel0>  fimction  program,  <north>,  <south>,  <top>  and  <bottom> 
matrices  were  passed  to  the  right  side  of  the  equation.  Their  values  were  calculated  by 
using  the  temperature  values  fi'om  the  previous  iterations  and  the  resulting  values  were 
added  to  the  values  of  <right>.  The  diagonal  matrix  was  solved  by  using  Gauss 
elimination  method  and  the  values  of  Tw  ,  Tp  and  was  foimd  for  the  present  iteration. 

In  the  <old_templ0>  function  program,  it  was  assumed  that  the  plane  was  moving 
with  the  arc  source  together  (Figure  A5).  Here,  due  to  the  non-uniform  grid  spacing,  it 
was  necessary  to  determine  the  new  position  and  the  temperature  value  of  each  grid  point 
by  using  second  degree  polynomial  enterpolation  (Figure  A6).  The  values  of  7^ ,  and 
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were  found  by  using  the  old  temperatures  from  the  previous  time  step.  These 
temperatures  was  used  to  determine  the  new  grid  temperature  with  the  second  degree 
polynomial  enterpolation. 


j 

i  \ 

•  t 

.7 

1 

»  A 

1 

•  1 

• 

1 


Figure  A5  Moving  plane 


(a) 


(b) 


Figure  A6  Grid  positioning  by  using  second  degree  polynomial 
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The  function  program  <save_alllO>  saves  all  3-D  temperature  data  at  each  time 
step  by  overwriting  onto.  The  other  function  program  <save_thelO>  saves  top  surface 
temperatures  at  each  time  step  into  a  different  file  name. 
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APPENDIX  B.  PROGRAM  CODES 


THE  MAIN  PROGRAM  (aathlO.m) 


%  Need  to  define  the  function  below  to  be  able  to  use  mcc 
%  The  symbol  "al"  itself  has  no  significance  -  aathlO  is  the 
%  name  of  the  main  program, 
function  [al] =aathlO () 

clear  %  clear  all  variables  -  not  present  in  "aathlOgoon.m" 

close  all  %  close  all  open  figure  windows 

%  save_me  is  a  counter  that  is  incremented  after  each  time  step 
%  and  is  used  to  decide  how  often  to  save  the  calc  temperature 
%  data. 
save_me=0 ; 

%  save_cotinter  is  a  counter  that  is  also  incremented  with  each 
%  time  step  and  is  appended  to  the  file  name  into  which  the  temp 
%  data  is  stored. 

save_counter=0;  %  not  present  in  "aathlOgoon.m” 

%  Tinf  is  the  ambient  temp  in  deg  C 
Tinf=27; 

%  TO  is  the  initial  temp  in  deg  C 
T0=27; 

%  sil  is  the  stefan-boltzmann  const  (SI  units) 
sil=5.67E-8; 

%  epsilon  is  the  emissivity  of  the  surface 
epsilon=0; 

%  x-vel  component  of  torch  speed  (m/s) 
vel_x=0 ; 

%  y-vel  component  of  torch  speed  (m/s) 
vel_y= 0.004; 

%  deltat_t  is  the  time  step  (in  secs) 
delta_t=0 .001; 

%  Distance  source  moves  in  the  x-direction  in  each  time  step  (in  m) 
deltax=vel_x*delta_t ; 


%  Distance  source  moves  in  the  y-direction  in  each  time  step  (in  m) 
deltay=vel__y*delta_t  ; 

%  q_up  is  the  heat  flux  distribution  on  the  upper  surface  (W/m^2) 
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%  q_up  is  the  heat  flux  distribution  on  the  upper  surface  (W/m^2) 
%  imposed  by  the  torch/arc  (implemented  in  matrix  form) 
q_up=0 ; 


%%  q  is  a  generic  variable  representing  heat  flux  that  may  be 
%%  present  at  any  of  the  other  surfaces  and  can  be  included 
%%  in  the  nodal  equations 

%  q_  is  constant  heat  flux  through  the  surfaces 

q_west=0; 

q_east=0; 

q_top=0; 

q_bottom=0; 

q_north=0; 

q_south=0; 

%  c  is  the  cp  [J/kg-K] 
c=509.3; 

%  ro  is  the  density  [kg/m'^S] 
ro=7854; 

I  Convection  heat  transfer  coefficients  from  the  faces  [W/m^2-K] 

he==0; 

hw=0; 

hn=0; 

hs=0; 

ht=0; 

hb=0; 

%  The  distances  between  the  grid  point  points 

%%  back  is  the  portion  of  the  moving  grid  behind  the  point 

%%  of  location  of  the  source  [in  m] 

back=.05; 

%  b  is  the  coefft  in  the  exponential  gridding  function  (e''(-b’^x)) 
b=. 13; 

%  number  is  the  no.  of  grid  points  in  the  "back”  region 
nuinber=20; 


%%  Variables  "b"  and  "number"  above  can  be  varied  to  control  the 
%%  features  of  the  variable  grid,  such  as  change  in  grid  spacing, 

%%  etc. 

%%  height  is  a  vector  that  holds  the  coords  of  the  grid  points  in  %% 
the  back-region.  Note  that  the  spacing  between  these  grid 
%%  points  is  the  diff.  between  consecutive  height  entries  and  is 
%%  stored  in  "int_back" 
height (1) =0; 
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%  Note  that  int_back(l)  >  int__back (number)  ,  i.e.  in  reverse  order 
for  kk=l: number 

height (kk+1) =back*  (1-exp (-b*kk) ) ; 
int_back (kk) =height (kk+1) -height (kk) ; 

end 


%%  multiply  is  a  scaling  factor (close  to  1)  that  scales  height  and 
%%  int_back  in  such  a  way  that  it  height (number )  =  back. 
multiply=back/ ( sum (int_back) -int_back(l) /2) 
int_back=multiply*int_back; 

%%  Similar  to  back,  below  are  the  variables  front  (front  region  of  %% 

%%  moving  grid  ahead  of  source),  width  (width  of  plate)  and 

%%  thickness.  All  distances  in  m. 

front=0.02; 

width=0 .09; 

thickness=0 . 04 ; 

%%  int_front  is  the  grid  spacing  in  the  front  region.  Initially 

%%  set  to  0.  Then  transfer  from  int-back  one  by  one  until  s\im  of 

%%  int-front  entries  exceeds  front. 

int__front=0; 

counter=0; 

while  ( sum ( int_f rent ) -int_f ront (length ( int_f ront ) ) /2 ) <=f ront 


%  int_front  below  continually  changes  size  as  the  spacings  are 
%  added. 

int_front= [int_front, int_back (length (int_back) -counter)  ]  ; 
counter=counter+l ; 

end 

%  Set  first  member  which  was  initially  0  to  null. 

%  Meaningful  spacing  starts  only  from  the  2nd  member  of  int_front 
int  front (!)=[] ; 


%%  The  grid  distance  matrix  in  the  y  direction  is  y_intl  starting 
%%  from  the  top  of  front  down  through  the  source  all  the  way  to 
%%  the  bottom  of  back. 

%  fliplr  is  to  ge  the  desired  ordering  from  front  to  back. 
y_intl= [fliplr (int_front ) , fliplr (int_back) ] ; 

%  int__side  holds  grid  spacing  values  over  the  width. 

%  Same  logic  as  for  int_front 

int_side=0; 

count er=0 ; 

while  (sum (int_side) -int^side (length (int_side) ) /2 ) <= (width/2 ) 
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int_side= [int_side^ int_back (length (int_back) -counter) ] ; 
counter=counter+l ; 
end 

int  side  (!)”[]  ; 


%  The  grid  distance  vector  in  the  x  dirextion  is  x_intl  (y_intl 
%  above. 

x_intl= [f liplr (int_side) , int^side]  ; 

%  The  grid  distance  vector  in  the  z  dirextion  is  z_intl  (y_intl 

%  above) 

z_intl=0; 

count er=0/ 

while  {suin(z_intl)  -z^intl  (length  ( z_intl)  )  /2 )  <=thickness 

z_intl= [z_intl, int_back (length (int_back) -counter) ] ; 
counter=counter+l ; 


end 

z_intl(l)  =  []  ; 


figure (1) 

subplot (3,1,1) 
plot (x_intl ) 
grid  on 


subplot  (3,1,2) 
plot (y_intl ) 
grid  on 

subplot (3,1,3) 
plot ( z_intl) 
grid  on 

%  row  is  the  number  of  y  nodes  (along  the  length) . 
row  =  length (y_intl ) “1 ; 

%  col  is  the  no.  of  x-nodes  (across  the  plate  width), 
col  =  length (x_intl) -1; 

%  surface  is  the  no.  of  z-nodes  (across  the  thickness), 
surface  =  length (z_intl)-l; 
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The  value  of  teta  determines  the  numerical  scheme  as  follows: 


%  for  Fully  Implicit  Method;  teta=l 
%  for  Crank-Nicholson  Method;  teta=0.5 
%  for  Explicit  Method;  teta=0 

teta=l; 


For  a  a*b*c  three  dimensional  matrix 


%  a:num  rows,  b:num  cols,  c:num 


% 

Q. 

Size 

of 

new - 

— > 

(a*b*c) *3 

*6 

% 

Size 

of 

top - 

— > 

(a*b*c) *1 

% 

% 

size 

of 

bottom 

— > 

{a*b*c) *1 

% 

% 

size 

of 

north  - 

— > 

(a*b*c) *1 

% 

% 

size 

of 

south  - 

— > 

(a*b*c) *1 

% 

% 

g. 

size 

of 

right  - 

— > 

{a*b*c) *1 

o 

% 

size 

of 

temp  — 

— > 

(a*b) *c 

surfaces . 

TP,  TE,  TW  coefficients  for  the 
whole  grids. 

TT  coefficients  for  the  whole 
grids . 

TB  coefficients  for  the 
whole  grids. 

TN  coefficients  for  the  whole 
grids . 

TS  coefficients  for  the  whole 
grids . 

The  values  of  the  righthandside 
of  the  matrix. 

The  temperatures  of  the  grids. 


%  new  is  the  coefft  matrix  for  all  the  grid  points  in  the  3-D 
%  volume. 


%  1st  col  holds  the  coeffts  of  the  West  (W)  nodes. 

%  2nd  col  holds  the  coeffts  of  the  (P)  nodes. 

%  3rd  col  holds  the  coeffts  of  the  East  (E)  nodes. 

new=zeros (row*col^surface,  3) ; 

%  right  is  the  const  (or  b-vector)  on  the  RHS 
right=zeros (row*col*surf ace, 1) ; 

%%  south,  north,  top,  bottom  hold  the  coeffts  of  the  corresponding  %% 

grid  points. 

south=right ; 

north=right ; 

top=right ; 

bottom-right ; 

%  Apply  the  initial  condition  to  the  plate. 

%  ”temp”  is  the  temperature  matrix. 

%  temp  is  a  matrix  that  holds  temperature  values  for  each  surface. 

%  Note  dimensions  of  temp  carefully. 

%  This  initialization  is  not  present  "aathlOgoon .m" . 
temp=TO*ones (row*surface, col) ; 
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%  Below  are  parameters  that  specify  the  heat /arc  source. 

%  rO  is  the  effective  radius  [m]  of  the  gaussian  heat  source. 
r0=0.0045; 

%  ddd  is  a  coefft  that  determines  the  gaussian  spread. 

%  The  larger  the  ddd  value,  the  more  pointed  the  source. 
ddd=3; 

%  qrate_total  is  the  total  energy  transfer  rate  from  the  source 
%  [Watts] . 
qrate_total  =  2544; 


%  qO  is  the  peak  heat  flux  of  the  gaussian  [W/m^2] 
q0=qrate__total*ddd/pi/r0'^2 ; 

%  Calculate  the  q_up  (heat  flux  distribution  on  the  top  face) 

%%  xcenter  and  ycienter  are  the  grid  element  numbers  of  the  grid 
%%  center  at  which  the  source/arc  is  located, 
xcenter^length (int_side) ; 
ycenter=length (int_front ) ; 


for  nn=l:row 
for  mm=l:col 

%  x(inm)  is  array  with  the  x-coords  of  the  mesh, 
if  mm  <  xcenter 

X (mm) =-sum (x_intl ( (mm+1 ) : xcenter) ) ; 
elseif  mm==xcenter 
X  (mm)  =0; 
else 

X (mm) =sum (x_intl ( (xcenter+1) :mm) ) ; 

end 

%  y(inm)  is  array  with  the  y-coords  of  the  mesh, 
if  nn  <  ycenter 

y (nn) =sum (y_intl ( (nn+1 ) : ycenter) ) ; 
elseif  nn==ycenter 
y (nn) =0; 
else 

y (nn) =-sum (y_intl ( (ycenter+1) :nn) ) ; 

end 

%  initialize  q_up  using  x()  and  y()  coord  info. 

q_up  (nn,mm)  =-q0*exp  (“ddd*  (x(mm)  '^2+y(nn)  ^2)  /r0^2)  ; 

end 

end 

figure 

%  a  mesh  command  to  view  the  generated  grid, 
mesh (x, y, q_up) 
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for  time=delta_t :  delta__t :  200 

%  200  is  the  final  time  in  secs  upto  which  the  program  should 

%  go .  .  . 

save__me-save_me+l; 

disp(*the  interpolation  time  is  *) 
tic  %  beginning  of  tic-toc  loop 

[temp] =old_templ0 (row, col, surface, temp, x, y, deltax, deltay) ; 

toe  %  end  of  tic-toc  loop. 


%  figure 

%  contour (x, y, temp (1 : row, 1 : col) ) 
disp(’The  matrix  formation  time  is  :’) 
tic  %  beginning  of  tic-toc  loop, 

%%  Start  to  form  the  matrix  by  defining  every  point. 

%  Beginning  of  creation  of  coefft  matrix 

%  loops  step  through  each  point  on  all  the  surfaces,  cols  and  rows 

for  kk=l: surface 
for  jj=l:row 
for  ii=l:col 


%  The  distances  between  the  reference  point  and  the  adjoining 
%  nodal  points. 

%  Distance  between  P&E  and  P&W 
xpe=x_intl (ii+1)  ; 
xwp=x_intl (ii)  ; 

%  Distance  between  P&N  and  P&S 
ypn=y_intl ( j j ) ; 
ysp=y_intl { j j+1)  ; 

%  Distance  between  P&T  and  P&B  ^ 

2pt=2_intl (kk) ; 

2bp=2_intl (kk+1)  ; 

%%  The  lateral  surface  areas  of  the  faces  of  the  control  voliome 
%%  around  the  reference  point  P. 

%  East-face  area. 

Ae= (ypn/2+ysp/2) * { zpt/2+2bp/2 ) ; 

Aw=Ae ; 
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%  North-face  area. 

An=(xpe/2+xwp/2) * ( 2pt/2+2bp/2 ) ; 

As=An; 

%  Top-face  area. 

At= (xpe/2+xwp/2) * (ypn/2+ysp/2 )  ; 

Ab-At ; 

%  The  si2e  of  the  control  volume. 

deltaV= (xpe/2+xwp/2) * {ypn/2+ysp/2) * ( 2pt /2+zbp/2 ) ; 


%  Routine  to  set  thermal  conductivities  below: 

%  Initiali2e  TEMP ( ) s  to  0  for  logical  ifs  below: 
TEMP(1)=0;TEMP(2)=0;TEMP(3)=0;TEMP(4)=0;TEMP(5)=0;TEMP(6)=0/ 


%  TPO  is  the  temperature  of  the  current  node. 

TP0=temp { (kk-l) *row+j j , ii) ; 

% 

% 

%  T_west=TEMP(l)  T_east=TEMP (2 )  T_north=TEMP (3) 

%  T_south=TEMP (4 )  T_top=TEMP ( 5 )  T_bottom=TEMP ( 6 ) 

% 

% 

if  ii==l  %  first  col,  no  west  neighbor,  so. . . 

TEMP(1)=TP0;  %  set  T_west  to  TPO 

elseif  ii==col  %  last  col,  no  east  neighbor,  so . 

TEMP(2)-TP0;  %  set  T_east  to  TPO 

end 

if  jj==l 

TEMP (3)=TP0; 
elseif  jj==row 
TEMP(4)=TP0; 

end 

if  kk=-l 

TEMP(5)=TP0; 
elseif  kk==surface 
TEMP(6)=TP0; 

end 

%%  If  current  node  P  has  a  valid  neighbor,  then  use  temperature 
%%  of  that  node .... 

if  (TEMP(1)==0) 

TEMP (1) =temp ( (kk-l) *row+j j  ,  ii-1)  ; 

end 
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if  (TEMP(2)==0) 

TEMP (2) =temp ( (kk-1) *row+j j , ii+1) ; 

end 

if  (TEMP(3)==0) 

TEMP  (3)  =teinp  (  (kk-1)  *row+j  j-1,  ii)  ; 
end 

if  (TEMP(4)==0) 

TEMP (4 ) =temp ( ( kk-1 ) *row+j j+1, ii) ; 
end 

if  (TEMP(5)==0) 

TEMP (5) =temp ( ( kk-2) *row+j j  ,  ii)  ; 
end 

if  (TEMP{6)==0) 

TEMP  (6)  =temp  (kk*row+j  j  ,  ii)  ; 
end 


%  Do  loop  to  calculate  k-values  of  all  grid  points. 

for  kjl=l : 6 

%  Now  set  the  thermal  conductivity  values. 
kelvin=TP0+273 . 15; 

kl=53;  %  nominal  kl-value  in  W/m-K  set  here,  but... 

%  Actual  kl-values  are  calc  below, 
if  (kelvin>=300  &  kelvin<1000) 
kl=53-0.04* {kelvin-300) ; 
elseif  (kelvin>=1000  &  kelvin<1800) 
kl=25+6.25E-3* (kelvin-1000) ; 
elseif  kelvin>=1800 
kl-125; 

end 

%  Find  the  thermal  conductivity  at  the  adjacent  point 
kelvin=TEMP ( kj 1 ) +273 . 15 ; 

k2=53;  %  nominal  k2-value  in  W/m-K  set  here,  but... 

%  Actual  k2-values  of  neighbors  calc  here, 
if  (kelvin>=300  &  kelvin<1000) 
k2=53-0.04* {kelvin-300)  ; 
elseif  (kelvin>=1000  &  kelvin<1800) 
k2=25+6.25E-3* (kelvin-1000)  ; 
elseif  kelvin>=1800 
k2=125; 

end 
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%  Harmonic  mean  thermal  conductivity  found  next 
NEW_K(kjl)  =  2*kl*k2/(kl+k2) ; 

end 

%  Now  redefine  k’s  of  the  neighboring  faces  with  values  from  NEW  K 

kw=NEW_K ( 1 ) ; 

ke=NEW_K(2) ; 

kn=NEW_K(3) ; 

ks=NEW_K ( 4 ) ; 

kt=NEW_K ( 5 ) ; 

kb-NEW_K(6) ; 

%  Coefficients  used  in  developing  the  discretised  equations. 

ae=ke*Ae/xpe; 

aw=kw*Aw/xwp; 

an=kn*An/ypn; 

as=ks*As/ysp; 

at=kt*At/2pt; 

ab=kb*Ab/zbp; 

apO=ro*c*deltaV/delta_t ; 

ap=teta* (ae+aw+an+as+at+ab) +apO; 

%  Coefficients  used  for  the  wall-medium  interface  part  of  the 

%  discretised  equations. 

coef fe=he*xpe/ (2*ke+he*xpe) ; 

coef fw=hw*xwp/ (2*kw+hw*xwp) ; 

coef fn=hn*ypn/ (2*kn+hn*ypn) ; 

coeffs=hs*ysp/  (2*ks+hs’^ysp)  ; 

coefft=ht*zpt/ (2*kt+ht*zpt) ; 

coef fb=hb*zbp/ (2*kb+hb*zbp) ; 

%  Coefficients  used  for  the  heat  flux  part  of  the  discretised 
%  equations. 

fluxcoef fe=xpe/ (2*ke+he*xpe) ; 
fluxcoef fw=xwp/ (2*kw+hw*xwp) ; 
fluxcoef fn=ypn/ (2*kn+hn^ypn) ; 
fluxcoef fs=ysp/ {2*ks+hs*ysp) ; 
fluxcoef ft=2pt/ (2*kt+ht*zpt ) ; 
fluxcoef fb=zbp/ {2*kb+hb*zbp) ; 

%  Radiation  coefft  for  surface  b.c.  -  see  thesis... 
radcoef f=sil*epsilon*2pt/ (2*kt+ht*zpt ) ; 


%  row_num  is  the  index  of  the  current  point  in  the  matrices  being  % 
used.  row_num  varies  from  1 : surface*col*row 
row_num= (kk-1 ) *row*col+ ( j j-1) *col+ii; 

%  For  lef t-back-top  corner 


if  ({kk==l)  &  (jj=-l)  &  (ii-=l)) 
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%  Temperature  values  for  current  time  step  for  each  relevant  node. 
TP0=temp ( (kk-l) *row+j j  ,  ii)  ; 

TE0=temp ( (kk-1) *row+j j  ,  ii+1)  ; 

TS0=temp ( (kk-1) *row+j j+1,  ii)  ; 

TB0=t emp ( kk*row+ j  j , ii ) ; 

%  Coeffts  of  grid  points;  note  south  and  bottom  later  are  moved  to 
%  the  rhs .  ■ 

%  "new"  is  the  A-matrix  in  the  solution. 

%  Note  teta  is  the  parameter  which  determines  the  numerical  method 
%  teta=l: fully  implicit;  teta=0 : explicit ;  teta=0 . 5 : Crank-Nicholson 
new(row_num, 2) =apO+teta* (ae+2*aw*coef fw+2*an*coef fn+as+2*at*coef ft+ab) 
new (row_num, 3) =-teta*ae; 
south (row^num, 1 ) =t eta* as; 
bottom (row_num, l)=teta*ab; 

%  "right"  is  the  constant  b-matrix  on  the  rhs... 
right {row_num, 1 ) = ( 1-teta) * (ae*TE0+as*TS0+ab*TB0 ) +  . . . 

(apO“ (1-teta) * (ae+2*aw*coef fw+2*an*coef fn+as+2*at*coef ft+ab) ) *TP0+  . . 
(2*aw*coef fw+2*an*coeffn+2*at*coefft ) *Tinf  - 
2*at*f luxcoef f t*q_up ( j j  ,  ii)  ... 

-2*aw*fluxcoef fw*q_west-2*an*fluxcoeffn*q_north-2*at*fluxcoef ft*q__top 

%%  radiation  &  q_boiling_coef f  belong  to  the  rhs  but  are  not  included 
%%  until  later  since  they  depend  on  the  current  temperature 


%  (a  form  of  quasi-linearization  of  nonlinear  temp  terms) 
radiation (row_num, l)=-2*at*radcoeff; 
q_boiling_coef f (row_num, 1) =-2*at*fiuxcoef ft ; 

%  for  right-back-top  corner 

elseif  ((kk==l)  &  (jj==l)  &  (ii==col) ) 

TP0=temp ( (kk-1) *row+j j  ,  ii) ; 

TW0=temp ( (kk-1) *row+j j , ii-1)  ; 

TSO^temp ( (kk-1) *row+j j+1,  ii)  ; 

TB0=temp  (kk^row-f  j  j  ,  ii)  ; 

new  (row__num,  2)  ==apO+teta*  (2*ae*coef fe+aw+2*an*coef fn+as+2*at*coef ft+ab) 
new  (row_num,  1)  =-teta*aw; 
south  (row_num^  1)  =t  eta*  as; 
bottom (row_num, 1 ) =teta*ab; 

right (row_num, l)=(l-teta) * (aw*TW0+as*TS0+ab*TB0) +  . . . 

(apO- ( 1-teta) * (2*ae*coef fe+aw+2*an*coef fn+as+2*at*coef ft+ab) ) *TP0+ 

(2*ae*coef fe+2*an*coef fn+2*at*coef ft ) *Tinf- 
2*at*fluxcoef ft*q_up ( j j  ,  ii) .  .  . 

-2*ae*fluxcoef fe*q_east-2*an*fluxcoef fn*q_north- 


2*at*fluxcoefft*q_top; 

radiation (row_num, 1) =-2*at*radcoef f  ; 
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q_boiling_coef f (row^num, 1 ) =-2*at ^fluxcoef ft ; 

%  for  left-front-top  corner 

elseif  ((kk==l)  &  (jj==row)  &  (ii==l)) 

TP0=temp ( (kk-1) *row+j j , ii) ; 

TE0=temp ( (kk-1) *row+j j,  ii+1)  ; 

TN0=temp ( (kk-1) *row+j j-1,  ii) ; 

TB0=temp (kk*row+j j / ii) ; 

new(row_num, 2)=ap0+teta* (ae+2*aw*coef fw+an+2*as*coef f s+2^at *coef ft+ab)  ; 
new(row_num, 3) =“teta*ae; 
north (row_num, 1 ) =teta*an; 
bottom (row_num, 1) =teta*ab; 

right (row_num, 1) = ( 1-teta) * (ae*TE0+an*TN0+ab*TB0) +  ... 

(apO- (1-teta) * (ae+2*aw*coef fw+an+2*as*coef f s+2*at*coef f t+ab) ) *TP0+  . . . 
(2*aw*coeffw+2*as’*^coeffs+2*at*coefft)  *Tinf-2*at*f  luxcoef  ft*q_up  ( j  j  ,  ii) 

-2*aw*fluxcoef fw*q_west-2*as*fluxcoeffs*q_south-2*at*f luxcoef ft^q  top; 
radiation (row_num, 1 ) =-2*at*radcoef f ;  ” 

q_boiling_coeff (row_num, 1) =-2*at*f luxcoef ft ; 


%  for  right-front-top  corner 

elseif  ((kk==l)  &  (jj==row)  &  (ii==col) ) 

TP0=temp ( (kk-1) *row+j j  ,  ii)  ; 

TW0=temp ( (kk-l) *row+j j,  ii-1)  ; 

TN0=temp  (  (kk-1)  *row-i- j  j-1,  ii)  ; 

TB0=t emp ( kk*row+ j  j , ii ) ; 


new  (row_num,  2)  =apO+teta*  (2*ae*coef  fe+aw+an+2*as*coef  f s4-2*at*coef  f t+ab)  ; 
new (row_num, 1 ) =“teta*aw; 
north {row_num, 1) =teta*an; 
bottom (row_num, 1) =teta*ab; 

right (row_num, 1) = (1-teta) * (aw*TW0+an*TN0+ab*TB0) +  ... 

(apO- (1-teta) * (2*ae*coef fe+aw+an+2*as*coef fs+2*at*coef ft+ab) ) *TP0+ 

(2*ae*coeffe  +2*as*coeffs  +2*at*coefft)*Tinf- 
2*at*fluxcoefft*q_up  { j  j  ,  ii)  ...  . 

-2*ae*f luxcoef f e*q_east-2*as*fluxcoef f s*q_south- 
2*at*fluxcoef ft *q_top; 
radiation (row_num,  1) =“2*at*radcoeff; 
q_boiling_coef f (row_num, 1) =-2*at*fluxcoef ft ; 
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%  for  left-back-bottom  corner 


elseif  ( (kk==surface)  &  {jj==l)  &  (ii==l)) 

TP0=temp ( (kk-1) *row+j j  ,  ii) ; 

TE0=temp ( (kk-1) *row+j j  ,  ii+1)  ; 

TS0=temp ( (kk-1) *row+j j  +  l, ii)  ; 

TT0=t emp ( ( kk-2 ) *row+ j  j  ,  ii )  ; 

new (row_num, 2) =apO+teta* (ae+2*aw*coef fw+2*an*coef fn+as+at+2*ab*coef fb) 
new  (row_niim,  3)  =-teta*ae; 
south (row_num, 1) =t eta* as; 
top (row_num, 1) =teta*at; 

right (row_num, 1 ) = ( 1-teta) * (ae*TE0+as*TS0+at*TT0) +  ... 

(apO- (1-teta) * (ae+2*aw*coef fw+2*an*coef fn+as+at+2*ab*coef fb) ) *TP0+  . . 
(2*aw*coef fw+2*an*coef fn+2*ab*coef fb) *Tinf  . . . 
-2*aw*fluxcoeffw*q_west-2*an*fluxcoeffn*q_north- 
2*ab*fluxcoef fb*q_bottom; 


%  for  right -back-bottom  corner 

elseif  ( (kk==surface)  &  (jj==l)  &  (ii==col) ) 

TP0=temp ( (kk-1) *row+j j  ,  ii)  ; 

TW0=temp ( (kk-1) *row+j j , ii-1) ; 

TS0=temp ( (kk-1) *row+j j+1, ii) ; 

TT0=temp( (kk-2) *row+j j,ii)  ; 

new (row_num, 2) =apO+teta* (2*ae*coef fe+aw+2*an*coef fn+as+at+2*ab*coef fb) 
new (row_num, 1) =-teta*aw; 
south  (row__num,  1)  =teta*as; 
top (row_num, 1) =teta*at ; 

right (row_num, 1) = (1-teta) * (aw*TW0+as*TS0+at*TT0) +  ... 

(apO- (1-teta) * (2*ae*coeffe+aw+2*an*coeffn+as+at+2*ab*coeffb) ) *TP0+ 

(2*ae*coeffe+2*an*coeffn+2*ab*coef fb) *Tinf  . . . 

-2*ae*fluxcoef fe*q_east-2*an*fluxcoef fn*q_north- 
2*ab*fluxcoef  fb*q_bottom; 


%  for  left-front-bottom  corner 

elseif  ( (kk==surface)  &  (jj==row)  &  (ii==l) ) 

TP0=temp ( (kk-1) *row+j j , ii) ; 

TE0=temp ( (kk-1) *row+jj  ,  ii+1)  ; 

TN0=temp  (  (kk-1)  *row+ j  j -1 ii)  ; 

TT0==temp  (  (kk-2)  *row+j  j  ,  ii)  ; 

new (row_num, 2) =apO+teta* (ae+2*aw*coef fw+an+2*as*coef fs+at+2*ab*coef fb) 
new (row  num, 3) =-teta*ae; 
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north (row_num, 1) =teta*an; 
top(row_num, 1) =teta*at; 

right  (row_nuin,  1)  =  (l-teta)  *  (ae*TE0+an*TN0+at^TT0 )  +  .  .  . 

(apO- (l”teta) * (ae+2*aw*coef fw+an+2*as*coef f s+at+2*ab*coef fb) ) *TP0+  . . 
(2*aw*coef fw+2*as*coef fs+2*ab*coef fb) *Tinf  .  .  . 

-2*aw*f luxcoef fw*q_west“2*as*fluxcoef fs*q_south- 
2’^ab*f  luxcoef  fb*q_bottom; 


%  for  right-front-bottom  corner 

elseif  ( (kk==surface)  &  (jj==row)  &  (ii==col)) 

TP0=temp ( (kk“l) *row+j j , ii) ; 

TW0=temp ( (kk-1) *row+j j, ii~l)  ; 

TN0=temp { (kk-l) *row+j j-1,  ii)  ; 

TT0=tenip  (  ( kk-2 )  *row+ j  j  ,ii); 

new(row_num, 2) -apO+teta* {2*ae*coef fe+aw+an+2*as*coef fs+at+2*ab*coef fb) 
new (row_num, 1 ) =-teta*aw; 
north {row^num, 1 ) =teta*an; 
top (row_num, 1 ) =teta*at ; 

right  (row_nuin,  1)  =  (1-teta)  *  {aw*TW0+an*TN0+at*TT0 )  +  ... 

(apO-  (1-teta)  *  (2*ae*coef fe4-aw+an+2*as*coef f s+at  +  2*ab*coef fb)  )  *TP0+ 

{2*ae’^coef  fe  +2*as*coeffs  +2*ab*coef  fb)  *Tinf  ... 

-2*ae*f luxcoef fe*q_east-2*as*f luxcoef fs*q_south- 
2*ab*fluxcoef  fb*q_bottoni; 


%  for  the  top  surface 

elseif  ((kk==l)  &  (jj>l)  &  (jj<row)  &  (ii>l)  &  (ii<col) ) 

TP0=temp ( ( kk-l ) *row+ j  j , ii ) ; 

TE0=temp ( (kk-l ) *row+j j , ii+1 ) ; 

TW0=temp ( (kk-l) *row+j j, ii-1) ; 

TN0=temp ( (kk-l) *row+j j-1, ii) ; 

TSp=temp ( (kk-l) *row+j j  +  1,  ii)  ; 

TB0==temp  (kk*row+j  j  ,  ii)  ; 

new (row^num, 2 ) =apO+teta* (ae+aw+an+as+2*at*coef ft+ab) ; 

new (row_num, 3) =-teta*ae/ 

new  (row^num,  1 )  ==-teta*aw; 

north  (row_nuin,  l)=teta*an; 

south  (row_num,  1 )  =teta*as; 

bottom (row^num, 1) =teta*ab; 

right (row__num, 1) = (1-teta) * (ae*TE0+aw*TW0+an*TN0+as*TS0+ab*TB0) +  ... 
(apO- (1-teta) * (ae+aw+an+as+2*at*coefft+ab) ) *TP0  +  2*at*coef ft*Tinf  - 

2*at*fluxcoef ft^q_up ( j j , ii) -2*at*fluxcoef ft*q_top; 
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radiation (row_num, 1 ) =-2*at*radcoef f ; 
q_boiling_coef f (row^num^ 1 ) =-2*at*fluxcoef ft; 


%  for  the  bottom  surface 

elseif  { (kk==surface)  &  (jj>l)  &  (jj<row)  &  (ii>l)  &  (ii<col) ) 

TP0=temp  {  (kk*“l)  *row+j  j  ,  ii)  ; 

TE0=temp ( (kk-1) *row+j j , ii+1) ; 

TW0=temp ( (kk-l) *row+j j , ii-1)  ; 

TN0=temp ( (kk-l) *row+j j-l,  ii) ; 

TS0=temp ( (kk-1) *row+j j  +  1,  ii) ; 

TTO^temp { (kk-2) *row+j j , ii)  ; 

new (row^num, 2 ) -apO+teta* (ae+aw+an+as+at+2*ab*coef fb)  ; 

new{row_num, 3) =-teta*ae; 

new  {row__num,  1)  =~teta*aw; 

north (row_num, 1) =teta*an; 

south (row_num, 1) -teta^as; 

top  (row___num,  1)  =teta*at; 

right (row_num, 1) = (1-teta) * {ae*TEO+aw*TWO+an*TNO+as*TSO+at*TTO ) +  . . . 
(apO“  (1-teta)  *  {ae+aw+an+as+at+2*ab*coef fb)  )  *TP0  +  2'^ab*coef fb*Tinf 

2*ab*f luxcoef fb*q_bottom; 


%  for  the  front  surface 

elseif  ((kk>l)  &  {kk<surface)  &  (jj==row)  &  {ii>l)  &  (ii<col) ) 

,TP0=temp ( (kk-1) *row+j j , ii) ; 

TE0=temp ( (kk-1) *row+j j , ii+1)  ; 

TW0=temp ( (kk-1) *row+j j , ii-1) ; 

TN0=temp ( (kk-1) *row+j j-1,  ii) ; 

TT0==temp(  (kk-2)  *row+j  j,ii)  ; 

•TB0=temp  (kk*row+j  j  ,  ii)  ; 

new(row_num, 2)=ap0+teta^ (ae+aw+an+2*as*coef f s+at+ab)  ; 

new (row_num, 3) =“teta*ae; 

new(row_num, l)=-teta*aw; 

north (row_num, 1) =teta*an; 

top (row_num, l)=teta*at; 

bottom (row_num, 1) =teta*ab; 

right (row_num^ 1 ) = ( 1-teta) * (ae*TEO+aw*TWO+an*TNO+at*TTO+ab*TBO) +  . . . 
(apO- (1-teta) * (ae+aw+an+2*as*coef f s+at+ab) ) *TP0  +  2*as*coef f s*Tinf 

2*as*fluxcoef fs*q_south; 


%  for  the  back  surface 

elseif  ((kk>l)  &  (kk<surface)  &  (jj==l)  &  (ii>l)  &  (ii<col)) 
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TPO=temp ( (kk"l) *row+j j  ,  ii)  ; 

TEO=temp ( (kk-1) *row+j j, ii+1) ; 

TWO=teinp  (  ( kk-1 )  ^row+ j  j  ,  ii-1 )  ; 

TSO=temp ( (kk-1) *row+j j  +  1,  ii) ; 

TTO^temp ( (kk-2) *row+j j  ,  ii) ; 

TBO==teinp  ( kk*row+ j  j  ,  ii )  ; 

new(row_num, 2) =apO+teta* (ae+aw+2*an*coef fn+as+at+ab) ; 

new(row_num, 3) =-teta*ae; 

new (row_num, 1 ) =-teta*aw; 

south (row_num, 1) =teta*as; 

top (row_num, 1) =teta*at; 

bottom (row^num, l)=teta*ab; 

right (row_num, 1) = (1-teta) * (ae*TE0+aw*TW0+as*TS0+at*TT0+ab*TB0) +  ... 
(apO- (1-teta) * (ae+aw+2*an*coef fn+as+at+ab) ) *TP0  +  2*an*coef fn*Tinf 

2*an*f luxcoef fn*q_north; 


%  for  the  left-lateral  surface 

elseif  ((kk>l)  &  (kk<surface)  &  (jj>l)  &  (jj<row)  &  (ii==l)) 

TP0=temp ( (kk-1) *row+j j , ii) ; 

TE0=temp ( (kk-1) *row+j j  ,  ii+1) ; 

TNO^temp ( (kk-1) *row+j j-1, ii) ; 

TS0=temp{ (kk-1) *row+j j+1, ii ) ; 

TT0=temp ( ( kk-2 ) *row+ j  j , ii ) ; 

TB0=temp (kk*row+j j , ii) ; 

new(row_num, 2) =apO+teta* (ae+2*aw*coef fw+an+as+at+ab) ; 

new (row^num, 3) =-teta*ae; 

north (row_num, 1) =teta*an; 

south (row^num, 1 ) =teta*as ; 

top ( row_num, 1 ) ~teta*at ; 

bottom (row_num, 1 ) =teta*ab; 

right {row_num, 1)= (1-teta) * (ae*TE0+an*TN0+as*TS0+at*TT0+ab*TB0) +  . . , 
(apO- (1-teta) * (ae+2*aw*coef fw+an+as+at+ab) ) *TP0  +  2*aw*coef fw*Tinf 

2*aw*f luxcoef fw*q_west ; 


%  for  the  right-lateral  surface 

elseif  ((kk>l}  &  (kk<surface)  &  (jj>l)  &  (jj<row)  &  (ii=-col)) 

TP0=temp { (kk-1) *row+j j  ,  ii)  ; 

TW0=temp ( (kk-l) *row+j j  ,  ii-1)  ; 

TN0=temp ( (kk-1 ) *row+j j-1,  ii)  ; 

TS0=temp ( (kk-1) *row+j j  +  1,  ii)  ; 

TT0=temp ( (kk-2) *row+j j  ,  ii) ; 

TB0=temp (kk*row+j j , ii) ; 
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new (row_num, 2 ) =apO+teta* (2*ae*coef fe+aw+an+as+at+ab) ; 

new  (row__num,  l)=“teta*aw; 

north  (row_nuin,  1 )  =teta*an; 

south ( row^num, 1 ) =teta*as ; 

top (row_num, l)=teta*at; 

bottom (row_num, l)=teta*ab; 

right (row_num, 1) = (1-teta) * (aw*TW0+an*TN0+as*TS0+at*TT0+ab*TB0 ) +  ... 
(apO- (1-teta) * (2*ae*coef fe+aw+an+as+at+ab) ) *TP0  +  2*ae*coef f e*Tinf 

2*ae*fluxcoef fe*q_east ; 


%  for  the  front-top  edge 

elseif  ((kk==l)  &  (jj==row)  &  (ii>l)  &  (ii<col)) 

TP0=temp { (kk-l) *row+j j , ii) ; 

TE0=temp ( (kk-l) *row+j j , ii+1) ; 

TW0=temp  (  (kk-l)  *row-i-j  j  ,  ii-1)  ; 

TN0=temp ( (kk-l) *row+j j-1, ii) ; 

TB0=temp (kk*row+j j , ii) ; 

new (row_num,  2) =apO+teta* (ae+aw+an+2*as*coef fs+2*at*coef ft+ab) ; 

new (row_num, 3) =-teta*ae; 

new {row_num, 1 ) =-teta*aw; 

north (row_num, l)=teta*an; 

bottom (row_num, l)=teta*ab; 

right (row_num^ 1) = (1-teta) * (ae*TE0+aw*TW0+an*TN0+ab*TB0) +  ... 
(apO- (1-teta) * (ae+aw+an+2*as*coef fs+2*at*coef ft+ab) ) *TP0  +  . . . 
(2*as*coeffs+2*at*coefft )  "^Tinf  -  2*as*f luxcoef fs*q_south- 
2*at*fluxcoef ft*q_top  ... 

-2*at*fluxcoef  ft*q__up  { j  j  ,  ii)  ; 
radiation  (row_num,  1)  =-2*at*radcoef  f ; 
q_boiling_coef f (row_num, 1 ) =-2*at*f luxcoef ft ; 


%  for  the  front -bottom  edge 

elseif  ( (kk==surface)  &  (jj==row)  &  (ii>l)  &  (ii<col) ) 

TP0=temp ( (kk-l) *row+j j , ii) ; 

TE0=temp ( (kk-l) *row+j j , ii+1) ; 

TW0=temp  (  (kk-l)  *row+j  j  ii-1)  ; 

TN0=temp ( (kk-l) *row+j j-1, ii) ; 

TT0=temp ( (kk-2) *row+j j , ii) ; 

new (row_num, 2) =apO+teta* (ae+aw+an+2*as*coef fs+at+2*ab*coef fb) ; 
new(row_n\am,  3)=-teta*ae; 
new  (row_num,  1)  =-teta*aw; 

north  (row_num,  1)  =teta*an; 
top (row_num^ 1) =teta*at ; 
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right  (row__num,  l)  =  (l-teta)  *  (ae*TEO+aw*TWO+an*TNO+at *TTO )  +  .  .  . 
(apO- (l~teta) * (ae+aw+an+2*as*coef f s+at+2*ab*coef fb) ) *TPO  +  . . 
(2*as*coeffs+2*ab*coeffb)  *Tinf  -  2’^as*fluxcoeffs*q_south- 
2*ab*fluxcoef fb*q_bottom; 


%  for  the  back-top  edge 

elseif  {(kk==l)  &  (jj==l)  &  (ii>l)  &  (ii<col)) 

TP0=temp ( ( kk“l ) *row+j  j  ,  ii ) ; 

TE0=temp ( (kk-1) *row+j j  ,  ii+1)  ; 

TW0=temp ( (kk-1) *row+j j , ii-1) ; 

TSO^temp ( (kk-1) *row+j j+1, ii) ; 

TB0=temp ( kk*row+ j  j,ii); 

new(row_num, 2) =apO+teta* (ae+aw+2*an*coef fn+as+2*at*coef ft+ab)  ; 

new (row__num,  3)  =-teta*ae; 

new (row_num, 1 ) =“teta*aw; 

south (row^num, 1) =teta*as; 

bottom (row_num, 1) =teta*ab; 

right  (row_num,  1)  =  (1-teta)  *  (ae*TE0+aw'^TW0+as*TS0+ab*TB0)  +  .  .  . 
(apO- (1-teta) * (ae+aw+2^an*coef fn+as+2*at *coef ft+ab) ) *TP0  +  . . 
(2*an*coef fn+2*at*coef ft ) *Tinf  -  2*an*fluxcoeffn*q_north- 
2*at*fluxcoefft*q_top . . . 

-2*at*f luxcoef ft*q_up ( j j , ii) ; 
radiation (row^num, 1) =-2*at*radcoef f ; 
q_boiling_coef  f  (row_num,  1)  =-2*at'^fluxcoef  f  t  ; 


%  for  the  back-bottom  edge 

elseif  ( (kk==surface)  &  (jj==l)  &  (ii>l)  &  (ii<col)) 

TP0=temp ( (kk-1) *row+j j,  ii) ; 

TE0=temp ( (kk-1 ) *row+j j  ,  ii  +  1 ) ; 

TW0=temp ( (kk-1) *row+j j  ,  ii-1)  ; 

TS0=temp ( (kk-1) *row+j j  +  1,  ii) ; 

TT0=temp ( ( kk-2 ) *row+j j  ,  ii) ; 

new(row_num, 2) =apO+teta* {ae+aw+2*an*coef fn+as+at+2*ab*coef fb)  ; 

new (row_num, 3) =-teta*ae; 

new(row_num, 1) =-teta*aw; 

south (row_num, 1) =teta*as; 

top (row_num, 1 ) =teta*at ; 

right  (row__num,  1)  =  (1-teta)  *  (ae*TEO+aw*^TWO+as^TSO  +  at*TTO)  +  ... 
(apO- (1-teta) * (ae+aw+2*an*coef fn+as+at+2^ab*coef fb) ) *TP0  +  . . 
(2*an*coef fn+2*ab*coef fb) *Tinf  -  2*an*fluxcoeffn*q_north- 
2*ab*f  luxcoef  fb'*'q_bottom; 
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%  for  the  front-left  edge 

elseif  ((kk>l)  &  (kk<surface)  &  (jj==row)  &  (ii==l) ) 

TP0=temp ( (kk-1) *row+j j , ii) ; 

TE0=terap  (  (kk-1)  *row+j  j  ,  ii+1)  ; 

TN0=temp ( (kk-1) *row+j j-1, ii) ; 

TT0=temp ( ( kk-2 ) *row+ j  j  ,  ii )  ; 

TB0=t emp ( kk*row+ j  j , ii ) ; 

new  (row__num,  2)  =apO+teta*  {ae+2*aw*coef fw+an+2’^as*coef fs+at+ab) 

new  (row_n\im,  3)  =-teta*ae; 

north  (row_num,  1 )  =teta*an; 

top (row_num, l)=teta*at; 

bottom (row_num, 1) -teta*ab; 

right {row_num, 1) = (l-teta) * {ae*TE0+an*TN0+at*TT0+ab*TB0) +  . . . 
(apO- (1-teta) * {ae+2*aw*coef fw+an+2*as*coef fs+at+ab) ) *TP0  +  . 
(2*aw*coef fw+2*as*coef fs) *Tinf  -  2*aw*fluxcoeffw*q_west- 
2*as*fluxcoef fs*q_south; 


%  for  the  front-right  edge 

elseif  ((kk>l)  &  (kk<surface)  &  (jj==row)  &  (ii==col) ) 

TP0=temp ( ( kk-1 ) *row+ j  j , ii ) ; 

TW0=temp { (kk-1) *row+j j , ii-1) ; 

TN0=temp ( (kk-1) *row+j j-1, ii) ; 

TT0=temp ( (kk-2) *row+j j  ,  ii) ; 

TB0=temp ( kk*row+ j j , ii ) ; 

new (row__num, 2) =apO+teta* (2*ae*coef fe+aw+an+2*as*coef fs+at+ab) 

new(row_mam, 1) =-teta*aw; 

north {row_num^ 1) =teta*an; 

top (row_num, 1) =teta*at ; 

bottom (row_num, l)=teta*ab; 

right (row_num, 1 ) = (1-teta) * (aw*TW0+an*TN0+at*TT0+ab*TB0) +  ... 
(apO-  (1-teta)  *  (2*ae’^coef  fetaw+an+2*as*coef fs+at+ab)  )  *TP0  +  . 
{2*ae*coef fe+2*as*coef f s )  *Tinf  -  2*ae*fluxcoeffe*q__east- 
2*as*f luxcoef fs*q_south; 


%  for  the  back-left  edge 

elseif  ( (kk>l)  &  (kk<surface)  &  (jj==l)  &  (ii==l) ) 

TP0=temp ( (kk-1) *row+j j , ii) ; 

TE0=temp ( (kk-1) *rowtj j , ii+1) ; 

TS0=temp ( (kk-1) *row+j j+1, ii) ; 

TT0=temp ( ( kk-2 ) *row+ j  j , ii ) ; 

TB0=temp ( kk*row+ j  j , ii ) ; 
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new (row_num, 2) =apO+teta* (ae+2*aw*coef fw+2*an*coef fn+as+at+ab) 

new  (row_nuiu,  3)  =-teta*ae; 

south (row_num, l)=teta*as; 

top(row_nuin,l)=teta*at; 

bottom (row_num, 1) =teta*ab; 

right (row_num, 1 ) = (1-teta) * (ae*TE0+as*TS0+at *TT0+ab^TB0 ) +  . . . 
(apO- (l“teta) * (ae+2*aw*coef fw+2*an*coef fn+as+at+ab) ) *TP0  +  . 
(2*aw*coef fw+2*an*coef fn) *Tinf  ~2*aw*fluxcoeffw*q_west- 
2*an*fluxcoef fn*q_north; 


%  for  the  back-right  edge 

elseif  {(kk>l)  &  {kk<surface)  &  {jj==l)  &  (ii==col) ) 

TP0=temp ( ( kk-1 ) *row+ j  j  ,  ii ) ; 

TW0=temp { ( kk-1 ) *row+j  j  ,  ii-1 )  ; 

TS0=temp ( (kk-1) *row+j jt 1,  ii) ; 

TT0=temp ( ( kk-2 ) *row+ j  j  ,  ii )  ; 

TB0=temp {kk*row+j j  ,  ii) ; 

new(row_num, 2) =apO+teta* (2*ae*coef fe+aw+2*an*coef fn+as+at+ab) 

new(row_num, 1) =-teta*aw; 

south (row_num, 1 ) =teta*as; 

top'(row_num,  1)  =teta*at  ; 

bottom (row_num, 1 ) =teta*ab; 

right (row_num, l)  =  (l-teta) * (aw*TW0+as*TS0+at*TT0  +  ab^TB0) +  .  .  . 
(apO- (1-teta) * (2*ae*coef fe+aw+2*an*coef fn+as+at+ab) ) *TP0  +  . 
(2*ae*coef fe+2*an*coef fn)  *Tinf  -  2*ae^fluxcoef fe*q__east- 
2*an*f  luxcoef  fn*q_north; 


%  for  the  left-lateral-top  edge 

elseif  ((kk==l)  &  (jj>l)  &  (jj<row)  &  (ii==l)) 

TP0=temp { (kk-1) *row+j j  ,  ii) ; 

TE0=temp ( (kk-1) *row+j j  ,  ii+1)  ; 

TN0=temp  (  (kk-1)  ’*’row+j  j-1,  ii)  ; 

TS0=temp ( (kk-1) *row+j j+1, ii) ; 

TB0=temp {kk*row+j j  ^ii); 

new (row_num, 2 ) =apO+teta* (ae+2*aw*coef fw+an+as+2*at*coef ft+ab) 

new (row_num, 3) =-teta*ae; 

north (row_num, 1 ) =teta*an; 

south (row^num, l)=teta*as; 

bottom (row_num,  l)==teta*ab; 

right (row_num, 1)= (1-teta) * (ae*TE0+an*TN0+as*TS0+ab*TB0) +  . . . 
(apO- (1-teta) * (ae+2*aw*coef fw+an+as+2*at *coef ft+ab) ) *TP0  +  . 
(2*aw*coef fw+2*at*coef ft ) *Tinf  -  2*aw*f luxcoef fw*q_west- 
2*at*fluxcoef ft*q_top. . . 

-2*at*f luxcoef ft *q_up ( j j , ii) ; 
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radiation  (row_num^  1)  =’'2*at*radcoeff ; 
q_boiling_coef f (row^num, 1) =-2*at*fluxcoef ft ; 


%  for  the  left“lateral~bottom  edge 

elseif  ( (kk==surface)  &  (jj>l)  &  (jj<row)  &  (ii==l)) 

TP0=teinp  (  (kk-1)  *row+j  j  ,  ii)  ; 

TE0=temp ( (kk-1 ) *row+j  j , ii+1)  ; 

TN0=teinp  (  (kk-1)  *row+j  j-1,  ii)  ; 

TS0=temp ( (kk-1) *row+j j+1, ii) ; 

TT0=temp ( (kk-2) *row+j j  ^ ii)  ; 

new  (row^num,  2 )  =apO+teta*  (ae+2*aw*coef fw+an+as+at+2*ab*coef fb)  ; 

new (row_num, 3) =-teta*ae; 

north (row_num, 1 ) =teta*an; 

south (row_num, 1) =teta*as; 

top  (row__niain,  1 )  =teta*at ; 

right (row_num, 1 ) = (1-teta) * (ae*TE0+an*TN0+as*TS0+at*TT0) +  ... 
(apO- (1-teta) * (ae+2*aw*coef fw+an+as+at+2*ab*coef fb) ) *TP0  +  . . 
(2*aw*coef fw+2*ab*coef fb) *Tinf  -  2*aw*fluxcoef fw*q_west- 
2*ab*fluxcoef  fb*q_bottom; 


%  for  the  right-lateral-top  edge 

elseif  ((kk”l)  &  (jj>l)  &  (jj<row)  &  (ii==col)  ) 

TP0=temp ( (kk-1) *row+j j , ii) ; 

TW0==teinp  (  (kk-1)  *row+j  j  ,  ii-1 )  ; 

TN0=temp ( (kk-1) *row+j j-1, ii)  ; 

TS0=temp ( (kk-1) *row+j j  +  1,  ii)  ; 

TB0=temp ( kk*row+ j  j , ii ) ; 

new(row_nuin,  2)  =apO+teta*  (2*ae*coef fe+aw+an+as+2*at*coef ft+ab)  ; 

new (row_num, l)=-teta*aw; 

north(row_nuin,l)=teta*an; 

south  (row__nuni,  1)  =teta*as; 

bottom (row_num, l)=teta*ab; 

right (row^num, 1)  =  (1-teta) * (aw*TW0+an*TN0+as*TS0+ab*TB0)  +  ... 
(apO- (1-teta) * (2*ae*coef fe+aw+an+as+2*at*coef ft+ab) ) *TP0  +  . . 
(2*ae*coef fe+2*at*coefft ) *Tinf  -  2*ae*fluxcoef fe*q_east- 
2*at*fluxcoef ft*q_top. . . 

-2*at*fluxcoef ft*q_up ( j j , ii) ; 
radiation (row_num, 1) =-2*at*radcoef f ; 
q_boiling_coef f  (row_num,  1 )  =-2*at'^fluxcoef  ft ; 


%  for  the  right-lateral-bottom  edge 

elseif  ( (kk==surface)  &  (jj>l)  &  (jj<row)  &  (ii==col) ) 
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TPO=temp  (  (kk-1)  *row+j  j,  ii)  ; 

TWO=temp  (  {  kk-1 )  *row+ j  j  ,  ii-1 )  ; 

TNO=temp  (  (kk-1)  *row+j  j-1,  ii)  ; 

TSO=temp  (  (kk-1)  *row+j  j  +  1,  ii)  ; 

TTO=temp  (  (kk-2)  *row+j  j,  ii)  ; 

new(row_num,  2)  =apO+teta*  ( 2*ae*coef f e+aw+an+as+at+2’^ab*coef fb)  ; 

new (row^num, 1 ) =-teta*aw; 

north (row_num, 1 ) =teta*an; 

south (row^num, 1) =teta*as; 

top (row_num, 1) =teta*at; 

right (row^num, l)=(l-teta) * (aw*TW0+an*TN0+as*TS0+at*TT0) +  ... 
(apO-  (1-teta)  *  (2*ae*coef  fe+aw-fan+as  +  at+2*ab*coef  fb)  )  *TP0  +  .  .  . 
(2*a€*coef fe+2*ab^coef fb)  *Tinf  -  2*ae’^f luxcoef fe*q_east- 
2*ab*f  luxcoef  fb*q_bottom; 


%  for  the  interior  nodes 

%%  as  expected,  note  that  the  interior  node  has  all  of  its  neighbors 

%%  contributing  to  the  coefft  matrices 

else 

TP0=temp ( ( kk-1 ) *row+ j j , ii ) ; 

TE0=temp ( (kk-1) *row+j j , ii+1) ; 

TW0=temp ( ( kk-1 ) *row+j  j , ii-1 ) ; 

TN0=temp ( (kk-1) *row+j j-1,  ii)  ; 

TS0=temp ( (kk-1) *row+j j  +  1,  ii)  ; 

TTO-temp ( (kk-2) *row+j j  ,  ii)  ; 

TB0=t emp ( kk*row+ j  j , ii ) ; 

new ( r ow_num, 2 ) =ap ; 
new (row__num,  3)  =-teta*ae; 
new (row^num, 1 ) =-teta*aw; 
north (row_num, 1 ) =teta*an; 
south (row^num, 1 ) =teta*as; 
top (row_num, 1 ) =teta*at; 
bottom (row_num, 1 ) =teta*ab; 

right (row_num, 1 )=( 1-teta) * (ae^TEO+aw^TWO+an^TNO+as^TSO+af'TTO+ab^TBO ) + 
(apO- (1-teta) * (ae+aw+an+as+at+ab) ) *TP0; 


end 

% 

the  "end*'  of 

if  loop 

to  decide  where  node  is  located 

end 

% 

the  "end"  of 

loop  for 

ii=  . 

.  .  (cols) 

end 

% 

the  "end"  of 

loop  for 

jj=  • 

.  .  (rows) 

end 

% 

the  "end"  of 

loop  for 

kk=  . 

.  .  (surfaces) 

toe 

% 

end  of  tic-toc 

loop 
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%%  solve  the  matrix  and  find  the  new  temperatures  by  using  the 
"aasolvelO .m” 

%%  function  -  row-by-row  sweeping  iterations  from  surface  to  surface 

[temp]  ==aasolvelO (row, col, surface, new, right, north, south, top, bottom, 

temp, radiation, q_boiling_coef f , Tinf , x_intl, y_intl, x, y) ; 

%  print  the  maximum  temp  in  the  entire  3-d  grid 
max  (max  (temp)  ) 


%%  save  the  row  number,  column_number,  surface  number,  time  and 
temperature 

%%  in  every  iteration 


save_me  %  print  value  of  save_me  counter 

%  time_step_interval  is  the  freq  at  which  to  save  temp  data 
time_step_interval=l  ; 

if  (save_me/time_step_interval)  ===  fix (save_me/time_step_interval) 

%%  save_counter  is  used  in  ”save_thelO .m”  function  to  apend  to  file¬ 
name 

%%  in  which  temp  data  is  being  saved 
save_counter=save_counter+l ; 

%  save_thelO  saves  top-surface  temp  at  each  time-step  into  a  diff  file 
name 

save_thelO (row, col, surface, time, temp (1 : row, 1 : col) , save_counter, b, x, y) ; 

%%  save_alllO  saves  all  3-d  temp  data  at  each  time-step  by  overwriting 
onto 

%%  the  same  file  name  -  imp  for  restarting  the  iterations... 

save_alllO (number, time, temp, b, delta_t , back, front , width, thickness, save_c 
ounter, save_me) ; 

end 


num_milli_secs==500;  %  freq  at  which  to  save  all  3-d  temp  data 
for  aa=l:20 

if  ( round  (save_me)  “round  (num_milli_secs’^aa)  ) 

save_thelO (row, col, surf ace, time, temp, save_counter, b, x, y) ; 

end 

end 

end 
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%  THE  FUNCTION  PROGRAM  (aasolvelO.m) 

%  this  program  solves  the  matrix  produced  by  the  main  program 
%  by  using  iterative  sweeping  method. 

% 

function  [temp]  =  aasolvelO (row, col, surface, new, right , north, south, top. 


bottom,  temp,  radiation,  q_boiling_coef  f ,  Tinf ,  x_intl,  y__intl,  x,  y)  ; 
%  for  a  a*b*c  three  dimensional  matrix. 


% 


size  of  new  - >  (a*b*c)*3 


size  of  top - > 

size  of  bottom  — > 

size  of  north  - > 

size  of  south  - > 

size  of  right  - > 


a*b*c) *1 
1 
1 
1 
1 


{a^b-^c)  * 
(a*b*c) * 
(a*b*c) * 
{a*b*c) * 


size  of  temp - >  (a*b)*c 

size  of  radiation  — >  (a*b)’ 
size  of  q_boiling_coef f  — > 
convection  heat  fluxes  from 


TP,  TE,  TW  coefficients  for  the 
whole  grids 

TT  coefficients  for  the  whole  grids 
TB  coefficients  for  the  whole  grids 
TN  coefficients  for  the  whole  grids 
TS  coefficients  for  the  whole  grids 
The  values  of  the  righthandside  of 
the  matrix 

The  temperatures  of  the  grids 
1  The  radiation  coefficients 
(a*b)*l  The  boiling  and  free 
the  top  surface 


%  initialize  the  temp_old,  temp_old  is  the  old  temperatures  of  the 
%  grid  points.  It  is  used  to  compare  the  temperatures  before  and  after 
%  the  iterations.  If  the  difference  is  less  than  0.1,  the  iteration 
%  is  enough. 


temp_old=temp+ 100* ones (size (temp) ) ; 
it_num=0; 

temp_saturation=100; 

delta_iteration=0 . 1 ; 
count_iterat ion=0 ; 

while  max (max (abs ( temp_old-temp) ) )  >  delta_iteration 
count_iteration  =  count_iteration+l ; 

if  count_iteration>delta_iteration*50+5 
delta_iteration=delta_iteration+0 . 1 
end 

max  (max  (abs  ( temp_old-temp)  )  ) 

it  num~it  num+1  %  iteration  number 
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disp  (max  (max  (temp)  )  ) 
t  emp__o  1  d= t  emp  ; 


%  The  variables  of  the  boiling  regimes 

g=9.81; 

csf-0.0132/ 

ne=l; 

hfg=2257E3; 
ro_l=957.85; 
ro__v=0 , 6; 

surface_tension=58 . 9E~‘3; 

Prandtl=1.75; 

miw_l=2 . 775E-4 ; 

miw_v=12 . 02E-6; 

cpl=4211; 

cpv=2029; 

epsi_s=0 . 82; 

k_liq=0.682; 

k_vap=0. 0249; 

alpha_l=0 . 169E-6; 

T_sat=100; 

T_liquid=27; 
sigma_s=5. 67E-8; 


pisi_pisi=2eros (row, col)  ; 


%  start  the  iterations 


for  cz=l: surface  %  for  the  surfaces- 

for  cy=l:row  %  for  the  rows 

indexl= (cz-1) *row*col+ (cy-1) *col+l; 
index2= (cz-1) *row*col+cy*col; 

newl=new ( indexl : index2 , : ) ; 
topl=top (indexl : index2, 1)  ; 
bottoml=bottom ( indexl : index2 , 1 ) ; 
northl=north ( indexl : index2 , 1 ) ; 
southl=south (indexl : index2, 1) ; 
rightl=right (indexl : index2, 1) ; 
radiationl==zeros  (length  (indexl :  index2) ,  1)  ; 
q_boiling_free=zeros  (length  (indexl :  index2)  ,1)  ; 

for  cx=l:col  %  for  the  columns 
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%  multiply  the  coefficients  by  the  temperatures 


if  cz==l 

qqqq=0;  %  heat  flux  to  the  outside  because  of  boiling  or  free 
convection 


%  The  radiation  heat  flux  from  the  top  surface 

radiationl ( 1 : col , 1 ) =radiation ( ( (cy- 
1 ) *col+l )  :  (cy*col ) , 1 ) . * (temp (cy, cx) ^4  )  ; 

radiationl (cx, 1 ) =radiation ( (cy-1) *ccl+cx, 1) * (temp(cy,cx) ^4) / 


%  The  boiling  or  free  convection  heat  flux  from  the  top  surface 

delta_temp  =  temp(cy,cx)  -  temp_saturation; 

if  (delta_temp<=5 ) 

%  Rayleigh  number 

T_film= (temp (cy, cx) +Tinf ) /2+273.15; 
g=9.81; 

beta=-0. 000000017937* (T_film^2) +0 . 0000188*T_film-0 . 0037515; 
alpha=-l  .  875E-12*  (T_film^2) +1.5525E-9*T_film-1.4995E-7; 
nu=l  .  125E-10*  (T_film"2)  -8 . 285E-8 *T_f  ilm+1  ..5584E-5; 
kconduc=-0. 000008125*  (T_f  ilm"^2  )  +0 . 0063775*T_f  ilm-0 . 568  95  ; 

%  characteristic  length  L 

%  The  free  convection  coefficient  will  be  found  for  a  1*1  m2  area 
L=l/4; 


Ra=abs (g*beta* (temp (cy, cx) -Tinf ) * (L^3) / (nu*alpha) ) ; 
if  Ra<16 

Nusselt=l ; 


else 


if  (Ra<-1E7) 

Nusselt=0 . 54 * (realpow (Ra,  0 . 25) ) ; 
else 

Nusselt=0 .15* (realpow (Ra,  (1/3) )) ; 
end 


end 

h_boiling=Nusselt*kconduc/L; 
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q_boiling_free (cx, 1) =q_boiling_coef f ( (cy-l) *col+cx, 1) * .  .  . 

h_boiling* (temp (cy, cx) -Tinf ) ; 

qqqq=h_boiling* (temp (cy, cx) -Tinf) ; 

end 


if  (delta_temp>5  &  delta_temp<=30) 

q_prime_s=miw_l*hfg*sqrt  (g*  (ro_l-ro_v)  /surface-tension)  *.  .  . 
(cpl^deltS— temp/csf /hfg/Prandtl^ne) ^3; 

tau=pi/3*sqrt (2*pi) *sqrt (surf ace-tension/g/ (rO-l-rO-V) ) * . .  . 
(realpow(abs  (rO— v'^2/surface_tension/g/  (rO— 1- 

rO-V) ) , 0.25) ) ; 

q— boiling^ f ree (cx, 1) =q— boiling— coeff ( (cy- 
1) *col+cx, 1) *q— prime— s* ... 

(1+ (2*k— liq* (T—Sat-T— liquid) / . . , 

sqrt (pi*alpha-l*tau) ) *24/ (pi*hfg*rO-V) * . . . 

realpow (abs (rO— v^2/surface— tension/g/ (rO—1- 

rO-V)),.25)); 

qqqq=q— prime— s* . , . 

(1+ (2*k— liq* (T—Sat-T— liquid) / . . . 

sqrt (pi*alpha-l*tau) ) *24/ (pi*hfg*rO-V) * . . . 

realpow (abs (rO— v^2/surface— tension/g/ (rO— 1- 

ro  v) ) , .25) ) ; 


end 

if  (delta— temp>30  &  delta-temp<=120) 

q-max=0 . 14  9*hfg*rO— v*realpow (abs (surface— tension*g* (rO— 1- 
ro-v) /rO-V^2) ,0.25) ; 

q— min=0 . 09*rO— v*hf g*realpow (abs (surf ace— tension*g* (rO— 1-rO— v)  / .  .  . 
(rO-l+rO-V)^2)  ,0.25) ; 

1 

log_q=logl0  {qjxiax/q_jci±n)  /loglO  (30/120)  *loglO  (delta-temp/120)  tloglO  (q-mi 
n)  ; 

q-boiling— free (cx, 1 ) =q— boiling— coeff ( (cy-1) *col+cx, 1) *10^ (log-q)  ; 

qqqq=10'^  (log-q)  ; 

end 

if  (delta— temp>120) 

lambda=2*pi*realpow(abs (surfacO-tension/g/ (rd-l-rO-V) ) ,  0.5)  ; 


97 


h_conduct=0.59*realpow(abs  (g*  (ro_lro__v)  *ro_v*k  vap^3^.  . 
(hfg+0 . 68*cpv*delta_temp) /  (lambda*miw_v^delta_temp) ) , . 25) ; 

h_radiate=sigma_s*epsi_s* ( ( (temp(cy,cx) ) ^4- (T_sat) ^4) /. . . 

(temp (cy,  cx) -T_sat ) )  ; 

h_boiling=h_conduct+0 . 75*h_radiate/ 

q_boiling_free (cx, 1 ) =q_boiling_coef f ( (cy-1 ) *ccl+cx,  1) *  . . . 

h_boiling* (temp (cy, cx) -Tinf ) ; 

qqqq=h_boiling^ (temp (cy, cx) -Tinf) ; 
end 

pisi_pisi (cy, cx) ^qqqq; 
end 


if  cz>l 

topi (cx, 1) =topl (cx, 1) *temp ( (c2-2) *row+cy,  cx)  ; 
end 


if  cz<surface 

bottoml  (cx,  1)  =bottoml  (cx,  1)  *temp  (cz*row-f-cy,  cx)  ; 
end 

if  cy>l 

northl (cx, l)=northl (cx, 1) *temp( (cz-1) *row+cy-l,  cx) ; 
end 

if  cy<row 

southl (cx, l)=southl (cx, 1) *temp ( (cz-l) *row+cy+l,  cx) ; 
end 

end 


%  indexl . . . index2  shows  the  rows  in  the  new  matrix  that  the 
%  calculations  are  done 


if  cz==l 


right 1 (1 : col, 1) =rightl (1 : col, 1 ) +radiationl (1 : col, 1 ) +q_boiling_f ree (1 : co 

1,1); 

end 

right 1 (l:col, l)=rightl (l:col, l)+northl (l:col,  1)  ... 

+  southl (1 ; col, 1 ) +topl ( 1 : col, 1) +bottoml ( 1 : col,  1 )  ; 
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% 


solve  the  diagonal  matrix  by  using  Gauss  elimination  method 


for  i=l: col-1; 

newl ( i , 3 ) =newl ( i , 3 ) /newl (1,2); 
right 1 (i, 1) =rightl (i, 1) /newl (i,  2) ; 
newl  (i,  2)  =1; 

newl (i+1, 2)=newl (i+1, 2) -newl (i+1, 1) *newl (i, 3) ; 
rightl (i+l, 1) =rightl (i+l, 1) -rightl (i, 1) *newl (i+1, 1) ; 
newl (i+1, 1) =0; 


end 

rightl (col, 1) =rightl (col, 1) /newl (col, 2} ; 
newl (col, 2) =1; 


for  i=col:-l:2 

rightl (i-1, 1) =rightl (i-1, 1) -rightl (i, 1) *newl (i-1, 3) ; 
newl (i-1, 3) =0; 

end 


%  The  calculation  is  over,  update  the  temperature  matrix  for  the  new 
values 

temp ( (cz-1 ) *row+cy, : ) =rightl (1 : col,  : )  * ; 

end 

end 

end 

%  figure 

%  mesh  (x,  y, pisi__pisi)  ; 

%  pisi_pisi 

%  title ( *heat  * ) 
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%  THE  FUNCTION  PROGRAM  (oldjemp  1 0.m) 

% 

%%  The  function  program  (olci_templO .m)  finds  out  the  new  positions  of 
%%  the  moving  grid  points  at  the  new  time  step  and  the  old  temperature 
%%  values  of  the  new  grid  points  by  using  second  degree  polynomial 
%%  enterpolation . 

function  [temp] =old_templO (row, col, surf ace, temp, x, y, deltax, deltay) ; 

%l  This  (for)  loop  finds  out  the  x-coordinates  of  the  new  area  where 
%%  the  grid  points  move  after  the  time  delta_t. 

for  forl=l:col 

new_pos_x (1, fori) =x  (1, fori ) +deltax; 


count  =  1; 


while  ((x(l, count)  <  new_pos_x (1, fori ) )  &  (count<col+10) ) 
if  count==col 

count=col+20;  %  There  was  ”break"  command  here  before. 

end 

count =count  + 1 ; 

end 

if  count>=col+20 
count=col; 

end 


if  (deltax>0) 

if  count>2  %  deltax>0  means  if  the  source  goes  to 

the  right 

count  =  count-2; 
else 

if  count  >  1 

count  =count-l; 

end 

end 

end 

if  (deltax<=0)  %  deltax>0  means  if  the  source  goes  to 

the  right 

if  count==col 

count=count-2 ; 
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else 

if  count>l 
count  =count“l; 

end 

end 


end 


pos_x(forl)  =  county- 

end 

%%  This  (for)  loop  finds  out  the  y-coordinates  of  the  new  area  where 
%%  the  grid  points  move  after  the  time  delta_t. 

forl=0; 

for  forl=l:l;row 


new_pos_y { fori ) =y ( 1 y  fori ) +deltay ; 
count  =  1; 


while  ((y(lycount)  >  new_pos_y ( 1, fori) )  &  (count<row+10) ) 

if  count=-row 
count=row+20; 

end 

count=count+l ; 

end 

if  count>=row+20 
count=row; 

end 

if  (deltay>0) 

if  count=“row 
count=count-2 ; 

else 

if  count>l 

count  =  count~l; 
end 

end 

end 

if  •(deltay<=0)  %  deltax>0  means  if  the  source  goes  to  the 

%  right 


%  deltax>0  means  if  the  source  goes  to 
%  the  right 
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if  count>2 
count=count-‘2  ; 
else 

if  count===2 
count  =count  ~ 1 ; 
end 

end 


end 


po's_y{forl)  =  count; 


end 

%%  Now,  we  know  the  coordinates  of  the  area  where  the  grid  points 
%%  move. We  have  to  use  second  degree  polynomial  interpolation  to  find 
%%  the  temperature  values  of  the  new  grid  points.  Here,  we  use  the 
%%  temperature  values  from  the  previous  time  step. 

for  f orl=l : surface 


for  for2=l:row 
for  for3=l:col 
%  temperatures 

yO=temp( (forl-1) *row+pos_y (for2) , for3)  ; 
yl=temp { (forl-l) *row+pos_y (for2) +1,  for3)  ; 
y2=temp ( (forl-1 ) *row+pos_y (for2) +2, for3)  ; 

%  positions 

xO=y (pos_y ( f or2 ) ) ; 
xl=y(pos_y (for2)+l) ; 
x2=y (pos_y (for2) +2)  ; 

%  our  second  degree  polynomial  is  y=a*x^2+b*x+c 
%  a,  b,  c  are  the  coefficients.  Now  find  these  coefficients... 
b=(yO*  {x2'^2-xl'^2) +yl*  (x0'^2-x2^2 ) +y2*  (xl^2-x0^2)  )  /  ... 

(xO*  (x2^2-xl^2)  +xl*  {x0'"2-x2^2)  +x2*  (xl^2-x0^2)  )  ; 

if  abs (xO) -=abs  (x2) 

a==(y0-y2+b‘^  (x2-x0)  )  /  {x0"2-x2"'2)  ; 
else 

a=  (yO-yl+b*  (xl-xO)  )  /  (x0^2-xl'"2)  ; 

end 

C=y0-a*x0^2-b*x0; 
temp_old( (forl- 

1)  *row+for2,  f or3 )  =a*new_pos_y  ( for2 )  ^2*fb*new_pos_y  (for2)  +c; 
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%  [pol_coef f ] =polyfit (y (pos_y (for2 ) : (pos_y ( for2 ) +2) ) ,  . . . 

%  temp ( ( ( (forl-1) *row+pos_y (for2) ) : ( (fori- 

1) *row+pos_y (for2)+2) ) , for3) ’ ,2) ; 

%  temp_old(  (forl- 

1)  *row+for2,  for3)  =polyval  (pol_coeff,new_pos_y  (for2)  )  ; 

end 

end 


end 


for  forl=l : surface 


for  for2=l:row 
for  for3=l:col 


yO=temp_old ( (forl-1 ) *row+for2, pos_x (for3) ) ; 
yl=temp_old { (forl-1) *row+for2, pos_x (for3) +1)  ; 
y2=temp_old(  (forl-1) *row+for2, pos_x (for3) +2) ; 

xO=x (pos_x ( f or3 ) ) ; 

2<l=x  (pos_x  ( f  or3 )  +1 )  ; 
x2=x  (pos__x  ( f  or3 ) +2 )  ; 

b=(yO*  (x2^2-xl^2)+yl*  (x0"^2-x2'^2 ) +y2*  (xl^2-x0^2)  )/  ... 

(xO*  (x2^2-xl^2)  +xl*  (x0'"2-x2^2)  +x2*  (xl^2-x0"^2)  )  ; 

if  abs  (xO) -'=abs  (x2) 

a- (y0-y2+b* (x2-x0) ) / (x0^2-x2^2)  ; 
else 

a= (yO-yl+b* (xl-xO) ) / (x0"2-xl^2)  ; 
end 

c=y0-a*x0^2-b*x0; 
temp ( (forl- 

1 ) *row+for2, for3) =a*new_pos_x (for3) ^2+b*new_pos_x (for3) tc; 

%  [pol_coef f  ]  =polyfit  (x  (pos_x  (for3)  :  (pos___x  (for3)  +2)  )  ,  ... 

%  temp_old( (forl- 

1) *row+for2,  ( (pos_x (for3) ) :  (pos_x (for3) +2) ) ) ,  2)  ; 

%  temp ( (forl-1) *row+for2, for3) =polyval (pol_coef f , new__pos_x (for3) ) ; 


end 

end 


end 
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THE  FUNCTION  PROGRAM  (save  alllO.m) 


%%  the  function  program  {save_alllO .m)  saves  all  3-d  temp  data  at  each 
%%  time-step  by  overwriting  onto 

function 

[]=save_alllO (number, time, temp,b, del ta_t , back, front, width, thickness, sav 
encounter, save_me) ; 

save  save  alllO 
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%  THE  FUNCTION  PROGRAM  (save_thelO .m) 

%  the  function  program  (save_thelO)  saves  top-surface  temp  at  each 
%  time-step  into  a  diff  file  name 

%%  (aathlO.m)  and  (aathlOgoon.m)  can  not  be  compiled  without  putting 
%%  "%"  symbol  in  front  of  save after  compiling  erase  "%"  symbol  and 
%%  save  the  program  again. 

function  []=save_thelO (row, col, surface, time, temp,  save_counter,b,x, y) ; 
save  ( [ 'temperaturelO_' , num2str (save_counter) ] )  ; 
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